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Chapter  1 
Introduction 


1.1  Introduction 

The  structures  of  today  are  rapidly  increasing  in  size  and  complexity  compared  to 
those  of  even  a  few  years  ago.  Additionally,  with  monetary  constraints  equating  to 
weight  constraints,  the  structures  of  today  are  also  becoming  lighter  than  ever  before. 
While  “large  flexible  structures”  is  typically  a  term  reserved  for  space  structures,  it 
can  also  include  modern  airframes  and  even  buildings.  The  common  denominator  for 
large  flexible  structures  is  their  inherently  low  damping  and  stiffness  characteristics, 
which  necessitate  the  use  of  active  control  systems  for  vibration  and  shape  control. 

In  the  past,  the  structure  and  its  control  system  have  been  designed  independently. 
Structural  design  and  optimization,  and  control  system  design  and  optimization,  have 
both  been  areas  cf  separate  research,  each  progressing  vigorously  along  its  own  path. 
However,  spurred  on  by  recent  proposals  of  new,  large,  highly  constrained  space 
structures,  the  question  has  arisen  as  to  whether  an  integrated  structural/control 
design  procedure  might  not  be  more  appropriate. 

The  main  aim  of  this  work  is  to  investigate  this  problem,  setting  forth  and  demon¬ 
strating  a  design  methodology  with  which  practical  integrated  design  problems  can 
be  solved.  The  remainder  of  this  chapter  is  devoted  to  a  historical  overview  of  the 
control/structure  interaction  (CSI)  problem,  with  respect  to  the  problem  formula¬ 
tion  and  the  numerical  solution  techniques  employed.  Chapter  2  presents  the  general 
mathematical  programming  formulation  of  the  design  problem,  then  considers  how 
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the  particular  controller  methodology  under  consideration  modifies  and  adds  to  this 
probem  formulation. 

In  Chapter  3,  the  specific  numerical  solution  technique  used  in  this  work  is  pre¬ 
sented.  An  algorithm  is  developed  to  solve  the  mathematical  programming  problem 
as  set  up  in  Chapter  2,  and  a  numerical  technique  for  the  special  case  when  full  state 
feedback  control  is  used  is  considered.  In  Chapter  4,  the  solution  algorithm  is  illus¬ 
trated  on  a  representative  structure  —  the  DRAPER  I  Tetrahedral  Truss  structure. 
Finally,  from  these  results  some  general  conclusions  about  this  problem  are  drawn, 
and  recommendations  for  future  topics  of  research  given,  in  Chapter  5. 

1.2  Historical  background 

The  first  papers  actively  investigating  the  question  of  simultaneous  structure  and 
control  design  began  appearing  in  the  literature  around  19S3  [Han,  Kom,  Hal-1].  Since 
then,  there  has  been  a  growing  interest  in  this  subject  from  other  authors,  although 
the  field  itself  is  still  in  relative  infancy.  This  section  will  attempt  to  review  the  current 
state  of  the  art  in  the  simultaneous  structure/control  optimization  of  large  flexible 
structures.  The  section  is  broken  into  two  main  parts  dealing  with  the  general  problem 
formulations  presented  and  the  numerical  solution  techniques  employed.  Although 
these  are  somewhat  related,  they  are  largely  separate  since  several  numerical  methods 
could  be  used  to  solve  each  formulation,  and  several  formulations  could  be  solved  using 
the  same  numerical  procedure. 

1.2.1  Problem  Formulation 

Most  authors  begin  by  assuming,  if  only  implicitly,  that  the  simultaneous  design  of 
the  structure  and  the  control  system  will  definitely  reap  benefits  that  outweigh  the 
increased  computational  expense  that  will  inevitably  follow.  On  the  surface,  this 
would  seem  a  reasonable  assumption,  since  a  traditional  sequential  design  process 
may  involve  several  iterations  by  both  the  structural  design  group  and  the  control 
synthesis  group.  However,  Haftka  et  al.  [Haf-1]  found  that,  specifically  in  his  case  for 
the  integrated  thermal/structural  design  problem,  integrating  normally  independent 
design  methodologies  does  not  always  justify  the  effort  and  expense  involved. 
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Therefore,  to  investigate  the  possible  benefits  of  an  integrated  structure/control 
design  approach,  Ilaftka  et  al.  [Haf-3]  consider  small  changes  in  the  structural  pa¬ 
rameters.  They  find  that  significant  enhancements  in  the  control  can  be  obtained  by 
these  small  changes.  That  is,  small  structural  modifications  can  result  in  large  control 
effort  reductions  to  meet  the  same  vibration  control  requirements.  This  result,  and 
other  recent  studies  also  demonstrating  synergistic  effects  of  an  integrated  design  (for 
example,  [Kho-l,Sal,Mes]),  are  seen  as  a  justification  for  further  investigation  into 
the  problem. 

Using  a  conventional  design  approach  for  a  controlled  structure,  one  would  first 
optimize  the  structure  alone,  then  design  a  control  system  for  this  baseline  structure. 
This  process  may  then  be  iterated  until  both  the  structure  and  control  system  meet 
necessary  constraints  and  objectives.  The  structure  may  have  constraints  such  as  a 
maximum  total  mass  and  minimum  open-loop  frequencies  and  spacings,  stress  and/or 
displacement  constraints  under  various  loading  conditions,  and  perhaps  even  speci¬ 
fied  open-loop  mode  shapes.  The  control  system  constraints  may  include  closed-loop 
damping  and  frequency,  frequency  spacing,  total  available  control  energy,  and  maxi¬ 
mum  output  response  constraints.  In  this  conventional  approach  then,  the  structure 
would  be  des’gned  to  minimize  the  total  mass  subject  to  open-loop  structural  con¬ 
straints,  and  the  controller  would  then  be  designed  in  some  optimal  fashion,  perhaps 
as  an  linear  quadratic  gaussian  (LQG)  regulator,  subject  to  the  closed-loop  con¬ 
straints.  The  objective  to  be  minimized  in  the  control  design  would  be  a  quadratic 
form  of  the  type 


xT  Ax  +  utBu\  (1.1) 

Therefore,  a  “classical”  approach  to  the  simultaneous  structure/control  optimiza¬ 
tion  would  be  to  simultaneously  minimize  the  weighted  sum  of  the  total  mass  and  the 
above  quadratic  form,  subject  to  all  of  the  structural  and  control  constraints.  That 
is,  define  the  objective  as 

J{p)  =  ai\V(p)  +  a2  J  |xT A(p)x  +  uT B(p)u}  dt  (1.2) 

There  are  a  number  of  authors  who  formulate  the  integrated  control/structure 


/  [xtAx  +  utBuj  dt 


or 
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design  optimization  using  this  “classical”  approach.  Salama,  Hamidi  and  Demsetz 
[Sal],  Mossac,  Turner  and  Soosaar  [Mes],  Miller  and  Shim  [Mil-1]  and  Onoda  and 
Ilaftka  [Quo]  are  members  of  this  group.  In  addition,  these  authors  simplify  the 
problem  by  considering  that,  for  any  fixed  parameter  vector  p,  the  optimal  constant 
gain  control  is  given  by  the  well  known  Linear  Quadratic  Regulator  (LQR)  solution. 
Then,  equation  (1.2)  reduces  to 

J(p)  =  a^p)  +  a2xT{p)P3S{p)x(p)  (1.3) 

where  Pss  is  the  steady-state  solution  to  the  appropriate  Riccati  matrix  equation. 
This  significantly  reduces  the  dimension  of  p  since  the  control  gains  are  no  longer 
included  in  the  optimization,  and  can  result  in  a  significant  computational  savings. 

Hale,  Lisowski  and  Dahl  [Hal- 2]  also  define  an  objective  function  of  the  form  of 
(1.3),  and  augment  it  with  functions  of  auxiliary  coordinates,  Lagrange  multipliers 
and  adjoint  displacements.  Necessary  conditions  for  an  optimal  control  at  fixed  p, 
along  with  necessary  conditions  for  an  optimal  p,  can  then  be  derived  by  the  calculus 
of  variations.  The  resulting  coupled  nonlinear  equations  must  be  solved  numerically 
due  to  their  general  complexity. 

As  noted  by  Hale  and  Lisowski  [Hal-3],  one  of  the  major  problems  in  assuming  a 
quadratic  controller  penalty  function  ol  the  form  of  (1.1)  is  the  difficulty  in  assigning 
a  meaningful  relative  penalty  on  the  deflection  vs.  the  penalty  on  the  control  energy. 
This  same  problem  exists  in  the  separate  optimal  control  problem,  but  is  exacerbated 
in  the  integrated  design  problem  since  the  structural  parameters  are  varied.  Addition¬ 
ally,  we  must  also  assign  a  relative  weight  to  the  structural  mass  part  of  the  combined 
objective  (1.2)  vs.  the  control  part  of  (1.2),  i.e.  we  must  assign  values  to  aj  and  a?  in 
equation  (1.2).  The  optimal  design  so  calculated  will  be  very  sensitive  to  the  values 
of  these  weighting  parameters,  but  present  methods  rely  more  on  artwork  to  choose 
these  weights  than  on  any  formalized  algorithm.  One  way  to  avoid  choosing  these 
weights  would  be  to  extremize  only  one  part  of  the  joint  objective  function  (1.2),  and 
include  the  effects  of  the  other  part  explicitly  into  the  constraints. 


In  [Hal-3],  the  authors  use  the  total  maneuver  control  cost  as  the  objective  to  be 
minimized,  subject  to  initial  and  final  displacement  and  velocity  constraints  (maneu¬ 
ver  constraints).  Khot,  Eastep  and  Venkayya  [Kho-1],  Khot  [Kho-2],  Navid  [Nav]  and 
Manning  and  Schmit  [Man]  all  consider  the  problem  of  integrated  control/structure 
design  as  a  mass  minimization  subject  to  constraints  on  the  closed-loop  system  re¬ 
sponse  (displacement  and/or  eigenvalue  constraints).  While  the  control  cost  can  easily 
be  a  function  of  structural  parameters  through  the  system  matrices,  it  is  usually  as¬ 
sumed  that  the  total  mass  is  not  a  function  of  controller  parameters  (i.e.  neglecting 
the  mass  of  the  actuators  and  their  energy  source).  Therefore,  when  considering  a 
minimum  mass  design,  it  is  the  interaction  with  the  constraints  on  the  closed-loop 
system  that  provides  the  simultaneous  design. 

Bodden  and  Junkins  [Bod]  seek  to  minimize  the  norm  of  a  vector  of  output  feed¬ 
back  gain  matrix  elements  to  “herd”  the  closed-loop  eigenvalues  into  an  acceptable 
region  of  the  complex  plane.  In  this  sense,  the  norm  of  this  vector  was  taken  as  a 
measure  of  the  control  effort.  In  [Jun],  Junkins  and  Rew  extend  the  ideas  presented 
in  [Bod]  to  cases  where  the  controller  is  designed  as  a  classical  LQG  regulator.  Here, 
the  elements  of  the  weight  matrices  in  the  quadratic  performance  index  are  included 
in  the  parameter  vector  p,  and  the  authors  show  that  not  all  weight  matrices  are 
created  equal. 

Haftka,  et  al.  [Haf-3,Haf-2]  choose  a  similar  objective  function,  in  this  case  the 
sum  of  the  gains  of  the  actuators.  Since  in  their  case  the  actuators  act  as  viscous 
dashpots  attached  to  a  flexible  bar,  the  sum  of  these  gains  gives  a  measure  of  the 
total  amount  of  damping  supplied  by  the  actuators  and  colocated  velocity  sensors. 
In  this  sense,  the  sum  of  the  gains  is  also  a  measure  of  control  effort. 

Rather  than  minimize  the  sum  of  the  structural  mass  and  a  control  performance 
term,  Lim  and  Junkins  [Lim-1]  instead  review  three  measures  of  closed-loop  stability 
robustness  which  could  be  extremized.  This  paper  is  essentially  a  condensed  version 
of  the  Ph.D.  dissertation  of  Lim  [Lim-2],  Large  flexible  structures  have  many  closely 
spaced,  lowly  damped,  low  frequency  modes,  and  since  present  finite  element  tech¬ 
niques  will  only  accurately  calculate  the  first  few  modes,  there  may  be  considerable 
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uncertainty  in  modal  frequencies  and  mode  shapes  within  the  bandwidth  of  the  con¬ 
troller.  Therefore,  the  closed-loop  system  should  be  as  robust  as  possible  to  variations 
in  these  parameters. 

The  first  robustness  measure  is  an  eigenvalue  sensitivity  norm,  which  is  the  max¬ 
imum  singular  value  (L2  norm)  of  a  matrix  product  representing  an  upper  bound 
on  the  square  of  a  weighted  eigenvalue  error  norm  for  normalized  perturbations  p. 
This  measure  would  be  minimized  to  give  the  structure  stability  robustness.  Another 
measure  is  the  well-  known  condition  number,  which  can  be  related  to  the  radius  of 
uncertainty  within  which  all  eigenvalues  are  perturbed  due  to  an  error  in  the  system 
matrix.  At  the  very  least,  we  would  like  the  closed-loop  system  to  remain  asymp¬ 
totically  stable  in  the  presence  of  parameter  uncertainties,  and  this  corresponds  to 
maximizing  the  stability  robustness.  It  is  shown  that  this  is  equivalent  to  minimizing 
the  condition  number  of  the  closed-loop  system  matrix. 

The  third  robustness  measure  is  due  to  Patel  and  Toda  [Pat],  and  is  related  to 
the  spectral  radius  of  the  solution  of  an  associated  Lyapunov  equation.  This  mea¬ 
sure  should  be  maximized.  The  major  problem  with  measures  of  robustness  is  their 
conservatism.  However,  the  authors  in  [Lim-I]  conclude  from  their  numerical  studies 
that  a  significant  improvement  in  the  robustness  measure  does  indeed  correspond  to 
a  significant  improvement  in  the  actual  stability  robustness  of  the  closed-loop  system. 

Adamian  [Ada]  chooses  another  measure  of  robustness,  and  formulates  the  prob¬ 
lem  as  a  minimization  of  a  composite  objective  function  which  is  a  linear  combination 
of  the  structural  weight  and  the  robustness  measure.  In  his  work,  the  sensitivity  of 
the  closed-loop  eigenvalues  with  respect  to  plant  uncertainties  is  taken  as  the  measure 
of  robustness.  Closed-loop  eigenvalue  constraints  are  also  included. 

Gustafson,  et  al.  [Gus]  set  the  problem  up  in  a  unique  fashion.  They  have  de¬ 
veloped  a  method  to  obtain  equivalent  continuum  models  for  truss  type  structures 
with  a  pattern  of  repeating  elements.  From  this  equivalent  continuum  model  (dis¬ 
crete  modelling  is  used  for  rigid  components  of  the  structure),  PDE’s  are  obtained 
and  solved  in  the  frequency  domain  so  that  explicit  dependence  on  structural  param- 
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eters  is  achieved.  The  solution  is  in  distributed  form,  which  thus  retains  all  modal 
information  available  in  the  PDE.  The  sets  of  structural  and  control  parameters  can 
be  combined  into  a  single  optimization  problem  -  the  simultaneous  control/structure 
design  problem. 

Hale  [Hal-4],  proposes  an  approach  to  the  simultaneous  synthesis  for  vibration  sup¬ 
pression  using  the  set  theoretic  analysis  method  for  linear  dynamic  systems.  First, 
external  disturbances  are  assumed  to  be  unknown  but  bounded  time-varying  pro¬ 
cesses.  Second,  inequality  constraints  on  outputs  (measurements)  and  inputs  (con¬ 
trols)  are  included  explicitly.  The  purpose  of  the  integrated  synthesis  is  to  maximise 
the  amplitude  bound  on  the  unknown  disturbances  while  explicitly  satisfying  input 
and  output  constraints. 

A  similar  approach  is  taken  by  Slater  [Sla],  in  that  no  explicit  relative  weight 
between  the  structural  parameters  and  the  control  terms  is  required.  Rather,  the 
control  is  viewed  as  complementing  the  effect  of  the  structure  in  reducing  deformations 
due  to  an  applied  stochastic  disturbance  model.  The  formulations  by  Manning  and 
Schmit  [Man]  and  Lust  [Lust]  are  similar  except  that  the  excitation  forces  are  assumed 
deterministic  and  known,  and  are  explicitly  harmonic  in  the  case  of  Lust.  Thus,  in 
all  members  of  this  last  group,  real  physical  constraints  drive  the  solution  rather  than 
the  vague  quadratic  penalties  in  a  cost  functional. 

One  thing  to  keep  in  mind  concerning  the  preceeding  discussion  is  that  the  choice 
of  objective  to  be  extremized  will  usually  be  highly  mission  dependent.  The  problem 
formulation  will  always  involve  some  element  of  engineering  judgement. 

1.2.2  Problem  Solution 

Once  the  problem  has  been  formulated  as  a  nonlinear  programming  problem,  their 
solution  presents  a  challenge.  There  is  a  critical  need  for  a  dependable  optimization 
tool  that  efficiently  handles  a  large  number  of  nonlinear  inequality  constraints  and 
a  very  high  dimensional  design  space.  In  the  literature,  many  numerical  solution 
approaches  are  presented,  almost  all  of  which  are  gradient  based.  All  solution  tech- 
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niques  are  iterative  due  to  the  nonlinearity  of  the  problems  to  be  solved,  and  may  only 
converge  to  local  minima  due  to  the  fact  that  the  objective  function  will  usually  not 
be  globally  convex.  In  fact,  this  problem  of  convergence  to  local  minima  is  inherent 
to  gradient  based  approaches,  and  generally  cannot  be  avoided. 

The  gradients  of  the  objective  and  constraint  functions  with  respect  to  the  design 
parameters  will  be  required.  For  highly  nonlinear  and  large  dimensionality  problems, 
numerically  evaluating  these  gradients  using  finite  differencing  can  be  computationally 
expensive,  and  the  gradients  may  include  a  significant  error  over  their  analytically 
evaluated  counterparts.  Such  errors  may  halt  the  solution  algorithm  prematurely. 
Therefore,  if  at  all  possible,  efficient  analytical  expressions  for  the  gradients  should 
be  derived.  Even  so,  convergence  rates  may  be  slow,  hence  numerical  solutions  for 
even  relatively  simple  structures  may  still  be  computationally  expensive. 

For  algebraic  optimization,  where  there  are  no  differential  type  constraints,  there 
are  many  commercially  available  nonlinear  programming  codes  that  can  be  used  to 
solve  general  nonlinear  constrained  optimization  problems.  Haftka,  et  al.  [Haf-3, 
IIaf-2]  used  NEWSUMT  [Miu],  which  employs  an  extended  interior  penalty  function 
formulation  using  Newton’s  method  with  approximate  second  derivatives  for  each 
constrained  minimization.  NEWSUMT-A  [Gra-1]  was  developed  from  NEWSUMT 
by  adding  constraint  approximations  and  a  move-limit  strategy,  and  is  used  by  Onoda 
and  Haftka  [Ono], 

The  VMCON  optimization  subroutine  [Cra],  which  is  based  on  Powell’s  algorithm 
for  nonlinear  constraints  that  uses  Lagrangian  functions,  was  used  by  Khot,  Eastep 
and  Venkayya  [Kho-1]  and  Khot  [Kho-2].  Gustafson,  et  al.  [Gus]  used  the  IMSL 
subroutine  ZXSSQ  [IMSL].  This  subroutine  finds  the  least-squares  solution  based  on 
a  modification  of  the  Levenberg-  Marquardt  algorithm  which  eliminates  the  need  for 
explicit  derivatives.  This  is  a  special  case  because  the  example  considered  in  [Gus]  had 
no  constraints.  ZXSSQ  can  not  be  used  for  the  solution  of  constrained  optimization 
problems. 

Manning  and  Schmit  [Man]  used  the  subroutine  CONMIN  [Van-1]  to  solve  an  ap- 
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proximate  problem  that  is  derived  from  the  nonlinear  optimization  problem.  CON- 
MIN  was  also  used  by  Lust  [Lust]  and  Navid  [Nav]  to  solve  their  sequences  of  approx¬ 
imate  problems.  A  mathematical  programming  program  that  includes  CONMIN,  as 
well  as  a  variety  of  other  optimizing  algorithm  options,  is  ADS  [Van-3].  This  was 
used  by  Adamian  [Ada]  to  directly  solve  his  optimization  problem.  Both  CONMIN 
and  ADS  are  public  domain  codes,  widely  used  because  of  their  generality  and  ease 
of  application. 

Salama,  Hamidi  and  Demsetz  [Sal]  propose  an  iterative  solution,  based  on  a  modi¬ 
fied  Newton-Raphson  scheme,  of  the  appropriate  Kuhn-Tucker  optimality  conditions. 
However,  this  solution  requires  the  calculation  of  the  second  derivatives  of  the  objec¬ 
tive  function  with  respect  to  the  design  variables.  For  high  order  systems  especially, 
evaluating  these  second  derivatives  becomes  very  costly.  In  order  to  avoid  this  ex¬ 
pense,  Miller  and  Shim  [Mil-1]  explore  the  use  of  the  method  of  steepest  descent 
to  minimize  the  objective.  If  no  constraints  are  active,  this  method  is  the  same  as 
the  gradient  projection  method  mentioned  above.  When  constraints  are  active,  the 
objective  is  augmented  by  these  active  constraints  using  Lagrangian  multipliers,  and 
this  augmented  objective  is  then  minimized  by  steepest  descent  as  before. 

Rather  than  solving  the  full  nonlinear  programming  problem,  some  authors  prefer 
to  work  with  linearized  versions  of  the  full  nonlinear  equations.  In  this  sort  of  iterative 
algorithm,  move  limits  must  be  placed  on  the  design  variables  in  order  that  the  local 
linearity  be  satisfied.  Bodden  and  Junkins  [Bod]  and  Junkins  and  Rew  [Jun]  linearize 
the  objective  and  constraint  functions  about  the  current  value  of  the  design  variables. 
The  initial  design  point  is  taken  as  the  output  of  a  typical  sequential  design  process. 
The  authors  adopt  a  “minimum  modification'’  strategy,  whereby  from  the  infinity  of 
possible  solutions  to  a  small  number  of  equations  (linearized  objective  and  constraints) 
in  a  large  number  of  unknowns  (elements  of  the  design  variable  vector),  the  solution  is 
chosen  that  minimizes  the  changes  in  the  design  variables  while  achieving  the  desired 
changes  in  the  objective  and  constraints. 

The  numerical  proceedure  used  here  is  equivalent  to  a  gradient  projection  con¬ 
strained  optimization  algorithm.  Convergence  of  the  iterative  algorithm  can  be  en- 
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hanccd  by  employing  a  continuation  proceedure,  whereby  a  new  parameter  is  intro¬ 
duced  to  allow  the  creation  of  “stepping  stone”  problems  that  are  arbitrarily  close  to 
the  converged  neighbouring  solutions  [Bod,Hor].  The  continuation  approach  allows 
one  to  almost  arbitrarily  satisfy  the  local  linearity  requirements  by  choosing  a  small 
enough  step  in  the  continuation  parameter.  It  also  allows,  for  example,  restrictive 
constraints  to  be  applied  gradually  from  less  restrictive  ones. 


Lim  and  Junkins  [Lim-1]  propose  an  approach  to  the  solution  of  nonlinear  pro¬ 
gramming  problems  involving  sequential  linear  programming.  Here,  the  nonlinear 
problem  is  transformed  into  the  standard  linear  programming  problem  by  lineariza¬ 
tion,  a  translational  transformation  of  the  design  parameters,  and  by  adding  slack 
and/or  surplus  variables.  Once  in  this  standard  linear  format,  the  problem  formula¬ 
tion  is  straightforward,  and  very  efficient  and  reliable  linear  programming  software  is 
available  for  its  solution. 

Hale  [Hal-4]  introduces  the  use  of  set-theoretic  methods  for  the  solution  of  in¬ 
tegrated  structure/control  synthesis  problems.  Due  to  the  author’s  unique  problem 
formulation  mentioned  earlier,  the  constraints  on  the  output  and  control  terms  are 
box-type  constraints,  which  can  be  approximated  (conservatively)  as  ellipsoidal  sets. 
Set-theoretic  methods  can  then  be  applied  to  maximize  the  magnitude  of  the  dis¬ 
turbances  such  that  the  ellipsoidal  set  constraints  are  not  violated.  However,  the 
author  acknowledges  that  the  computational  burden  of  a  numerical  solution  using 
this  method  is  very  high.  This  would  prohibit  solving  problems  of  actual  structures, 
with  a  large  number  of  states  and/or  design  parameters,  using  this  method. 


For  problems  with  differential  or  integral  constraints,  a  calculus  of  variations  ap¬ 
proach  to  the  problem  solution  yields,  through  the  Fundamental  Theorem,  the  neces¬ 
sary  conditions  to  be  satisfied  by  an  extremal  [I\ir-1].  These  necessary  conditions  are 
the  Euler  equations,  and  in  general  are  nonlinear  two-  point  boundary  value  differen¬ 
tial  equations.  Generally,  these  problems  will  be  much  more  difficult  to  solve  than  a 
corresponding  algebraic  optimization  problem.  In  fact,  solving  these  equations  in  an 
iterative  framework  can  become  prohibitively  expensive  [Hal-2]. 
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One  method  of  solving  the  two-point  boundary  value  problems,  that  uses  a  pro¬ 
jected  gradient  iterative  method,  is  given  in  [Hal-2].  Similar  techniques  are  also  used 
by  Messac,  Turner  and  Soosaar  [Mes]  and  by  Hale  and  Lisowski  [Hal-3]  for  solution 
of  their  integrated  design  problem  formulation.  Basically,  the  iterations  proceed  in 
the  negative  direction  of  the  gradient  of  the  objective  (direction  of  greatest  decrease 
of  objective)  until  a  constraint  is  met.  At  this  point,  the  gradient  vector  is  projected 
onto  the  constraint,  and  this  becomes  the  direction  of  motion  for  that  iteration.  The 
comment  is  made  however,  that  these  algorithms  are  not  yet  efficient  enough  to  tackle 
large  problems  [Hal-3]. 

1.3  Scope  of  the  Present  Work 

The  work  contained  in  this  report  is  intended  to  extend  the  work  of  Slater  [Sla],  in 
which  a  preliminary  study  of  control/structure  optimization  was  performed.  Here, 
the  optimization  will  be  based  on  the  dynamic  response  of  a  structure  to  an  external 
disturbance  environment,  and  as  such  will  differ  from  most  previous  work.  Such  a 
response  to  excitation  approach  is  common  to  both  the  structural  and  control  design 
phases,  and  hence  represents  a  more  natural  control/structure  optimization  strategy 
than  relying  on  the  artificial  and  vague  control  penalties  used  by  other  authors.  The 
disturbances  are  to  be  considered  unknown  and  stochastic,  and  can  therefore  model  a 
wide  variety  of  actual  disturbance  states.  The  design  objective  is  to  find  the  structure 
and  controller  of  minimum  mass  such  that  all  the  prescribed  constraints  are  satisfied. 
The  controller  interaction  will  be  inserted  by  the  imposition  of  appropriate  closed-loop 
constraints,  such  as  closed-loop  output  response  and  control  effort  constraints. 
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Chapter  2 

The  Integrated  Control/Structure 
Optimization  Problem 


2.1  Introduction 

In  this  chapter,  the  integrated  control/structure  design  optimization  problem  is  for¬ 
mulated  as  a  general  nonlinear  constrained  mathematical  programming  problem.  This 
problem  formulation  is  quite  general  with  respect  to  the  constraints  that  can  be  ap¬ 
plied  and  the  controller  methodology  that  can  be  employed.  In  this  work,  only  the 
special  case  of  full  state  feedback  control  will  be  considered.  This  enables  the  problem 
to  be  partially  solved  analytically,  and  substantially  reduces  the  dimensionality  of  the 
numerical  problem  to  be  solved.  In  future  work,  more  realistic  and/or  sophisticated 
controller  methodologies  will  be  considered. 

2.2  General  Problem  Formulation 

The  integrated  control/structure  design  optimization  problem  can  be  stated  as  fol¬ 
lows:  find  the  vector  of  structural  and  controller  parameters  that  minimizes  the  mass 
of  the  structure  subject  to  limitations  on  the  available  control  energy,  and  a  set  of 
allowable  output  responses  to  a  set  of  prescribed  stochastic  disturbances.  Mathe¬ 
matically,  this  can  be  written  in  the  form  of  a  nonlinear  mathematical  programming 
problem  as 
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Minimize,  with  respect  to  p,  the  weight 
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is  the  jth  element  of  g,  an  m- vector  of  structural  constraints, 

is  an  n-vector  of  state  variables, 

is  an  nu- vector  of  control  forces, 

is  an  nw-vector  of  stochastic  disturbances, 

is  the  (n  x  n)  matrix  containing  the  system  dynamics, 

is  the  (n  x  nu)  matrix  containing  information  on  the 

locations  and  orientations  of  the  actuators, 

is  the  (n  x  nw)  matrix  containing  information  on  the  points 
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forces  involved  in  gCCi, 

Pi  is  the  maximum  allowable  value  of  the  ilh  expected  control 

effort  function 

R,  is  an  (nUi  x  nUi)  control  force  weighting  matrix, 

goc,  is  the  ith  output  response  constraint  cost  function, 

a,  is  the  maximum  allowable  value  of  the  ith  expected  output 

response  function  E[j/Jt 

Wi  is  an  (nd,  x  nd,)  output  response  weighting  matrix, 

yd>  is  the  ith  design  output  ndi -vector, 

Hd,  is  an  (rid,  x  n)  matrix  giving  the  relationship  between 

the  state  variables  and  ydi, 

Pi  is  an  iV- vector  of  minimum  design  variable  values,  and 

pu  is  an  N-vector  of  maximum  design  variable  values. 

The  side  constraints  are  the  strict  bounds  pt  and  pu  on  the  design  variables,  and 
are  vector  inequalities  that  are  imposed  element  by  element.  These  design  variable 
bounds  are  not  included  explicitly  as  constraints  in  the  problem  formulation.  Note 
that  the  structural  weight  and  the  structural  constraints  g  will  in  general  not  be 
functions  of  the  controller  design  variables  unless  the  controller  mass  is  included  in 
the  design.  Note  also  that  gcc,  is  a  weighted  mean  square  control  effort,  and  goc,  is 
a  weighted  mean  square  output  response.  Multiple  output  response  constraints  are 
allowed,  although  only  one  of  these  will  in  general  be  active  at  the  optimum  design. 
However,  all  of  the  control  effort  constraints  will  generally  be  active  at  the  optimum 
design. 

The  disturbances  acting  on  the  structure,  represented  by  w,  can  conceivably  be 
in  any  form  one  desires.  However,  it  is  assumed  in  this  work  that  w  is  a  stochastic 
disturbance,  specifically  zero  mean  Gaussian  white  noise  with  covariance  Xw.  The 
structure  will  respond  to  this  disturbance  with  some  transient  behaviour,  in  addition 
to  a  steady-state  response.  It  seems  reasonable  to  optimize  the  structure  for  the 
steady-state  response  to  the  disturbance  rather  than  the  transient  response  because 
the  transient  behaviour  will  normally  be  of  secondary  importance  to  the  response  ob- 
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jcctives  (such  as  long  term  pointing  accuracy).  Also,  for  steady  state  optimization,  the 
differential  equation  constraint  (state  equation)  can  be  replaced  with  a  steady  state 
covariance  equation,  and  the  mean  square  control  effort  and  output  response  con¬ 
straints  may  be  recast  in  terms  of  this  covariance.  Therefore  the  two-point  boundary 
value  problem  is  eliminated  and  the  numerical  solution  of  the  problem  is  significantly 
simplified. 

The  particular  form  of  this  covariance  equation  will  depend  on  the  form  of  the 
controller  used,  and  will  be  derived  for  each  controller  methodology  as  they  are  pre¬ 
sented.  Each  covariance  equation  will  be  of  the  general  form 

Fctx  +  XF5  +  Q  =  0  (2.2) 

where  Fci  is  the  closed-loop  dynamical  system  matrix,  X  is  the  symmetric  positive 
definite  covariance  matrix,  and  Q  is  a  symmetric  positive  definite  or  semi-definite 
matrix. 

To  obtain  the  first-order  necessary  conditions  for  optimality,  first  form  the  aug¬ 
mented  objective  function  Ja  as 

Ja  =  J(p)  +  tv[AI(FclX  +  XF];  +  Q)}+:jriK]gCCj(p)  + 

:= i 

+  £  K3oc,  ( P )  +  J2  mAp)  (2.3) 

;= l  J=1 

Then  the  Kuhn-Tucker  necessary  conditions  for  a  constrained  minimum  state  that, 
at  the  optimum  design,  the  following  conditions  must  be  satisfied  [Van-2]: 

>0  if  pi  =  pe, 

=  0  if  p(<  <Pi  <  pu, 

<0  if  p,  =  Pu, 
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dJg(p) 

dpi 
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0  for  j  = 

P:9j{p)  = 

0  for  j  = 

9 i(p)  <  0 
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9oCj  <  0 

for  j  -  1,.. 

•  •  ,n0 

9cc,  <  0 

for 

..,np 

FclX  +  XF?l+Q  =  0 
P'<P<Pu 


In  equations  (2.4)  -  (  2.8)  above,  Ax  is  the  matrix  of  Lagrange  multipliers  corre¬ 
sponding  to  the  covariance  equality  constraint  equation  (2.2),  and  A„,  Ay  and  pi  are 
the  vectors  of  the  Lagrange  multipliers  corresponding  to  the  inequality  constraints  of 
equations  (2.2L  Since  the  covariance  Lyapunov  equation  (2.2)  is  symmetric,  we  can 
assume  that  Ar  is  symmetric  also. 

Equations  (2.4)  indicate  that,  at  the  optimum  design,  the  gradient  of  the  objective 
function  augmented  with  the  constraints  via  Lagrange  multipliers  must  vanish  for 
those  variables  not  at  their  side  constraints.  This  derivative  can  be  written  a0 


OJaip) 

Op, 


+  JLwjF'.x  +  XFI  +  9)1  +  g 


+ 


Ogocjjp) 

dpi 


m 


+  Yip] 


OgAp) 

dpi 


for  i  =  1 , . . . ,  N 


(2.S) 
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These  gradients  are  taken  with  respect  tc  all  structural  and  controller  design 
parameters.  Equations  (2.6)  indicate  that  the  Lagrange  multipliers  must  be  non¬ 
negative  at  the  optimum  design.  For  those  inequality  constraints  not  active  at  the 
optimum,  equations  (2.7)  indicate  that  the  corresponding  Lagrange  multipliers  must 
be  zero.  This  is  because  non-active  inequality  constraints  should  take  no  part  in  the 
final  design.  Finally,  and  obviously,  the  optimum  design  should  satisfy  all  constraints, 
as  indicated  by  equations  (2.8). 

Now,  consider  how  the  general  formulation  as  presented  is  affected  by  the  partic¬ 
ular  choice  of  controller  methodology.  The  controller  constraints  take  on  particular 
functional  forms,  as  does  the  covariance  equation. 

2.3  Full  State  Feedback  Control 

The  simplest  form  of  feedback  control  is  to  feedback  the  entire  state  vector,  as  in 

u  =  -Kx  (2.9) 

where  K  is  the  (nu  x  n )  state  feedback  gain  matrix.  The  controller  design  variables 
for  this  case  will  be  the  nun  elements  of  K.  A  block  diagram  giving  the  structure  of 
this  controller  is  given  in  Figure  2.1.  Substituting  this  control  into  the  state  equation, 
and  assuming  that  the  disturbance  w  is  zero  mean  Gaussian  white  noise,  the  state 
covariance  matrix  X  for  this  case  can  be  found  from  the  Lyapunov  equation 

FclX  +  XFl  +  G.VXWGTW  =  0  (2.10) 

where  Fci  =  (F  —  GK)  is  the  stable  closed-loop  dynamical  matrix  for  the  full  state 
feedback  case,  X  —  E[xxr]  is  the  (n  x  n)  symmetric  state  covariance  matrix,  and 
Xw  =  EfuiuJ7’]  is  the  (nw  x  nw )  symmetric  covariance  matrix  for  the  stochastic  dis¬ 
turbances. 

Expressions  for  the  controller  constraints  in  terms  of  this  covariance  matrix  can 
then  be  obtained  as 
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(2.11) 


Qcc, 


tr[AT  RiKiX]  ,  ,  .  , 

=  - ai - 1  for  *  =  l.-.-.n/j 


ft 


9oc, 


ft 


for  i  =  1, . . .  ,n0 


(2.12) 


where  A',  is  the  (nUi  x  n)  partition  of  K  corresponding  to  U{.  It  is  assumed  that  the 
u,  are  independent,  and  that  u  is  ordered  as 


Ul 

w2 


(2.13) 


Note  that  nu,  =  nu,  and  that  the  columns  of  G  can  be  interchanged  to  force 
condition  (2.13)  to  be  satisfied.  The  elements  of  the  feedback  gain  matrix  K,  along 
with  the  elements  of  the  covariance  matrix  X,  are  the  independent  controller  design 
variables  in  the  augmented  Lagrangian  approach.  The  introduction  of  the  Lagrange 
multipliers  means  that  all  variables  axe  to  be  considered  independent,  even  though  X 
appears  dependent  on  a  particular  choice  of  K  through  equation  (2.10). 


Consider  equations  (2.S)  when  the  gradients  are  taken  with  respect  to  the  elements 
of  the  gain  matrix  K .  To  simplify  notation,  consider  a  matrix  A  with  elements  a^, 
and  define  for  scalar  s  the  matrix  by 


ds  _  ds 
dA  ..  da.ij 

L  J  t;  U 


(2.14) 


Using  this  notation,  the  Kuhn-Tucker  equations  for  differentials  with  respect  to  K 
can  be  written 
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dl< 


t  [tr  (A x[(F  -  GI<)X  +  X(F  -  GA')t])]  + 
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d 

dl< 


i=i  \ 


ft 


=  0  (2.15) 
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Let  R  be  aa  (nu  x  nu)  matrix  defined  in  a  block  diagonal  sense  as 


R  =  diag 


{GrW 


(2.16) 


then  equation  (2.15)  can  be  rewritten  as  (since  Au  is  not  a  function  of  K) 


0 

dl< 


-  [tr  (A*[(F  -  GI<)X  +  X(F  -  GI<f])  +  tr[I<TRKX}]  =  0  (2.17) 


Using  the  differentiation  rules  for  the  trace  operator  given  in  Appendix  A,  equa¬ 
tion  (2.17)  becomes 


-  2GrAtX  +  2RKX  =  0 


(2.18) 


which,  since  X  is  non-singular,  and  if  AUi  ^  0  for  all  i,  reduces  to 


K  =  R-'Gt\x 


(2.19) 


When  the  gradients  are  taken  with  respect  to  the  elements  of  X,  equations  (2.8) 
become 


d 

dX 


tr  (A*[(F  -  GI<)X  +  X(F  -  GK)t\ )  +  £  Au, 


i=i 


/tr  [KjRjKjX) 

\  ft 


-  1 


d 

dX 


EAv. 


Li=l 


^t  r{HjW,HdtX]  y 


=  0 


(2.20) 


Equation  (2.20),  along  with  the  definition  of  I\  in  equation  (2.19)  and  that  of  R  in 
equation  (2.16),  gives 


FT\X  +  A XF  -  AXGR "l GT Ax  +  W  =  0,  (2.21 ) 
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where 


(—I 

1=1  \ai  / 

which  is  an  algebraic  Riccati  equation.  Equations  (2.19)  and  (2.21)  define  the  solution 
to  the  optimal  control  problem 

min  Jc  =  J  |it1Vi  +  utRu  j  dt  (2.23) 

where  I\  is  the  optimal  steady-state  gain  matrix,  and  \x  is  the  steady-state  solu¬ 
tion  to  the  associated  Riccati  equation.  Sufficient  conditions  for  the  existance  of  a 
steady-state  Riccati  equation  solution  and  stability  of  the  control  law  given  by  equa¬ 
tion  (2.19),  is  that  the  system  be  completely  controllable  and  detectable.  Detectability 
indicates  that  all  unstable  modes  are  observable.  If  these  mildly  restrictive  conditions 
are  met,  the  closed-loop  system  Fci  is  guaranteed  to  be  asymptotically  stable. 

Therefore,  at  the  optimum  design,  where  the  Kuhn-Tucker  optimality  necessary 
conditions  for  a  minimum  will  be  satisfied,  the  control  design  variables  will  be  the 
solution  of  the  linear  quadratic  regulator  (LQR)  problem  equation  (2.23).  Although 
this  LQR  property  only  holds  true  at  the  optimum  point,  it  will  be  assumed  that 
at  every  point  in  the  design  cycle,  the  control  design  variables  will  be  found  as  the 
solution  to  the  optimal  control  problem  (2.23).  Therefore,  the  numerical  optimization 
problem  can  be  reduced  to  optimization  over  just  the  structural  design  variables,  along 
with  an  optimal  control  problem  solution  which  will  be  a  function  of  the  Lagrange 
multipliers  Au  and  Ay.  The  immediate  benefit  of  this  is  a  reduced  dimensionality 
nonlinear  programming  problem.  In  addition  since  under  the  conditions  stated  above, 
the  regulator  solutions  always  give  a  stable  closed-loop  system,  no  explicit  check  must 
be  performed  on  the  system  stability  during  the  solution  procedure. 

In  order  to  considerably  improve  the  conditioning  of  the  optimization  problem, 
normalize  the  objective  function  by  its  initial  value  —  the  initial  weight  of  the  struc¬ 
ture  W0  [Van-2].  If  the  structural  constraints  as  given  in  g  are  normalized,  and 
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since  the  controller  constraint  functions  have  been  written  in  normalized  form,  the 
optimization  objective  and  constraint  functions  will  then  all  be  of  order  one. 


Finally  then  for  this  case,  if  p4  is  the  vector  partition  of  p  representing  the  struc¬ 
tural  design  variables  only,  the  integrated  full  state  feedback  control /structure  design 
optimization  problem  can  be  stated  as  follows: 


Minimize,  with  respect  to  p4,  the  normalized  weight 


subject  to: 
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Figure  2.1:  Full  State  Feedback  Block  Diagram 


Chapter  3 

Numerical  Solution  Technique 


3.1  Introduction 


In  the  general  problem  formulation  presented  in  Chapter  2,  the  constraint  functions 
are  generally  highly  nonlinear  implicit  functions  of  the  design  variables.  Solution  of 
this  problem  could  be  attempted  by  the  direct  application  of  nonlinear  programming 
techniques;  that  is,  using  the  exact  functional  expressions  for  the  constraints.  How¬ 
ever,  this  approach  quickly  becomes  computationally  very  expensive  as  the  dimen¬ 
sionality  increases  since  the  full  objective  and  constraint  functions  must  be  evaluated 
at  every  step,  and  their  respective  gradients  at  most,  if  not  all,  steps  thoughout  the 
design  procedure.  Such  evaluations  tend  to  be  computationally  very  expensive. 

Approximation  techniques,  where  the  implicit  nonlinear  problem  is  replaced  by 
a  sequence  of  explicit  approximate  (although  not  necessarily  linear)  problems,  have 
been  shown  to  yield  efficient  and  powerful  algorithms  for  structural  design  optimiza¬ 
tion  (see,  for  example,  [Sch,Gra-2]).  It  has  yet  to  be  shown  however,  that  the  se¬ 
quential  approximations  technique  can  be  applied  successfully  to  the  integrated  con¬ 
trol/structure  design  optimization  problem.  Difficulties  arise  because  in  general,  the 
control  constraint  functions  tend  to  be  much  more  complex,  implicit  nonlinear  func¬ 
tions  of  the  design  variables  than  do  the  purely  structural  constraint  functions.  This 
means  that  to  use  approximation  techniques  in  the  integrated  control/structure  op¬ 
timization,  either  much  smaller  steps  in  the  design  space  will  be  required,  or  a  more 
sophisticated  approximation  scheme  must  be  implemented. 
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3.2  Solution  Algorithm 


In  this  method,  the  fully  constrained  nonlinear  optimization  problem  is  solved  by  the 
iterative  construction  and  numerical  solution  of  a  sequence  of  explicit  approximate 
problems.  The  approximate  problems  are  basically  first-order  Taylor’s  series  expan¬ 
sions  (in  some  convenient  intermediate  variables),  of  the  objective  and  constraint 
functions.  The  numerical  solution  is  accomplished  using  the  method  of  feasible  di¬ 
rections  as  implemented  in  the  mathematical  programming  code  ADS  [Kir- 1] .  We 
still  solve  the  approximate  problems  using  a  nonlinear  programming  code  because, 
depending  on  the  intermediate  variables  chosen,  the  approximate  functions  may  still 
be  nonlinear  functions  of  the  design  variables. 

A  schematic  of  the  solution  algorithm  is  shown  in  Figure  3.1.  The  solution  process 
begins  with  some  initial  structure,  which  is  analyzed  using  the  finite  element  tech¬ 
nique.  At  this  point,  the  gradients  of  the  active  constraint  set  are  evaluated,  and  the 
approximate  problem  is  formed,  with  respect  to  the  current  design.  Expressions  for 
the  gradients  of  all  constraints  considered  can  be  evaluated  analytically.  The  gradi¬ 
ents  for  constraints  common  to  all  controller  types,  along  with  those  introduced  for 
specific  controller  methodologies,  are  derived  in  a  later  section. 

The  approximate  problem  is  solved  with  ADS  using  an  active  constraint  set  strat¬ 
egy  to  reduce  the  dimensionality  of  the  approximate  problem  by  deleting  the  inactive 
constraints.  Move-limits  on  the  design  variables  are  imposed  during  the  solution  to 
ensure  that  the  design  remains  within  the  region  for  which  the  approximation  func¬ 
tions  are  of  acceptable  quality.  The  choice  of  move-limits  and  how  they  change  can 
have  a  significant  effect  on  convergence,  and  will  often  be  determined  from  numerical 
experience  with  the  particular  problem  at  hand. 

After  the  solution  of  the  approximate  problem,  the  structure  and  its  control  system 
are  deemed  optimal  if  a  convergence  test  on  either  the  absolute  or  relative  objective 
function  change  over  a  specified  number  of  successive  global  iterations  is  satisfied. 
Otherwise,  the  objective  and  active  constraint  gradients  are  evaluated  for  the  new 
design,  a  new  approximate  problem  formed,  and  the  process  above  is  repeated  in  an 


24 


iterative  manner.  The  solution  procedure  ends  when  the  design  variables  converge, 
or  when  the  number  of  iterations  exceeds  some  preset  maximum. 

Scaling  the  structure  and  controller  to  the  closest  constraint  surface  is  possible  in 
the  case  discussed  in  this  work,  because  of  the  special  assumed  form  of  the  controller. 
This  scaling  is  performed  on  the  initial  system,  and  on  the  system  immediately  fol¬ 
lowing  the  solution  to  the  approximate  problem.  The  scaling  procedure  is  discussed 
in  a  later  section.  The  solution  procedure  outlined  above  has  been  implemented  in 
a  research  computer  program  CSOPT,  written  in  FORTRAN,  and  run  on  a  CRAY- 
XMP  computer.  The  following  sections  outline  each  of  the  major  components  of  the 
solution  algorithm  in  more  detail. 


3.3  Scaling  to  the  Structural  Constraints 

The  structured  constraints  will  generally  be  confined  to  constraints  on  static  member 
stresses  or  nodal  displacements  for  specified  loading  conditions,  and/or  to  constraints 
on  the  open-loop  structural  frequencies.  Scaling  the  entire  structure  to  any  one  of 
these  constraints  is  generally  rather  straightforward  for  the  finite  elements  considered 
in  this  work. 

The  design  variables  for  a  truss/membrane/shear  element  type  structure  are  the 
elemental  cross-sectional  areas/ thicknesses.  Then  for  a  statically  determinate  struc¬ 
ture,  the  design  variables  are  inversely  proportional  to  the  static  displacements  and 
stresses.  For  indeterminate  structures,  this  inverse  proportionality  is  not  exact,  but  is 
still  a  good  approximation.  In  these  cases,  it  will  be  necessary  to  perform  the  inverse 
scaling  in  an  iterative  procedure,  but  only  a  few  iterations  may  be  required  to  achieve 
the  exact  scaling  desired,  depending  on  the  accuracy  required. 

Equal  scaling  of  all  the  structural  design  variables  will  only  affect  the  open-loop 
frequencies  if  the  structure  includes  some  nonstructural  mass.  Consider  the  Rayleigh 
quotient  expression  for  an  open-loop  frequency  A,  =  u>?  for  a  conservative,  non- 
gyroscopic  system,  as 
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(3.1) 


A  4>?K.4>j 
'  <f>J[Ms  + 

where  is  the  mode  shape  associated  with  At ,  and  Ms  and  M(  are  the  mass  matrices 
associated  with  the  structural  and  nonstructural  mass  respectively.  If  the  stiffness 
and  structural  mass  matrices  are  linear  functions  of  the  design  variables  (true  for 
truss,  membrane  and  shear  elements),  then  to  obtain  the  desired  open-loop  frequency 
Af,  consider  scaling  all  design  variables  by  the  same  scale  factor  s^,  to  get 

_  si4>Tj  Ks<t>j 
*'  ~  tflsiM.  +  Mtfc 

Note  that  the  nonstructural  mass  matrix  is  not  scaled.  If  the  ratio  of  actual 
open-loop  frequency  to  desired  open-loop  frequency  is 

h,  =  Y4  <3'3> 


(3.2) 


the  scaling  S{  necessary  to  obtain  this  ratio  can  be  obtained  by  substituting  equa¬ 
tions  (3.1)  and  (3.2)  into  equation  (3.3),  and  solving  for  the  scaling  parameter  s,,  to 
get 


_  _ _ 

'  4>J  [h,(Ms  +  Me)  —  M3]<f>{ 


(3.4) 


One  scaling  parameter  s,  for  each  frequency  constraint  will  be  obtained,  and  the 
scaling  parameter  of  the  most  violated  constraint  would  be  used  to  scale  the  system. 


3.4  Forming  the  Approximate  Problem 

Since  the  objective  is  linear  in  the  design  variables  (when  restricted  to  truss,  mem¬ 
brane  and/or  shear  finite  elements),  the  objective  gradient  is  constant  over  the  entire 
design  process  and  the  objective  function  can  be  evaluated  exactly  for  any  inter¬ 
mediate  design.  Therefore  the  objective  function  need  not  be  approximated  at  all. 
Approximation  functions  will  be  formed  for  both  the  structural  and  control  constraint 
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functions.  The  two  most  common  types  of  approximation  in  use  are  the  inverse  de¬ 
sign  variable  approximation  (see  [Can]  for  example)  and  the  hybrid  design  variable 
approximation  [Sta],  which  are  both  first-order  Taylor’s  series  expansions  in  some 
convenient  intermediate  design  variable. 

3.4.1  Inverse  Design  Variable  Approximation 

In  this  approximation  the  constraint  function  is  expanded  as  a  first-order  Taylor’s 
series  expansion  in  variables  that  are  the  inverse  of  the  design  variables.  As  noted 
in  a  previous  section,  this  approximation  is  exact  for  static  structural  constraints  on 
a  determinate  structure,  and  this  is  the  motivation  for  such  a  choice  of  intermediate 
variables  for  the  structural  design  variables.  Numerical  experience  gained  in  the 
current  work  has  shown  that  this  approximation  also  works  quite  well  on  the  controller 
constraints  in  some  cases. 

A  generic  constraint  approximation  function  written  in  terms  of  the  intermediate 
design  variables  is  then 


9i{y)  =  9i(Vo)  +  £  JT 
f=i  dVi 


(yj-Vjo) 


where  the  y}  =  are  the  (intermediate)  inverse  design  variables.  The  constraint  gra¬ 
dients  are  obtained  with  respect  to  the  direct  design  variables  p,  hence  equation  (3.5), 
written  in  terms  of  the  pj ,  becomes 


*M«*w-fy*f*  (£-£) 


Note  that  these  approximation  functions  are  nonlinear  in  the  design  variables  pj. 


The  gradients  of  these  approximation  functions  can  be  found  by  differentiating 
equation  (3.6)  with  respect  to  the  direct  design  variables  pj,  to  get 


Note  that  this  gradient  is  not  constant,  but  depends  on  the  current  design  variable. 


3.4.2  Hybrid  Design  Variable  Approximation 

The  use  of  approximation  functions  is  motivated  by  the  desire  to  obtain  a  low-order 
approximation  that  best  approximates  the  actual  constraint  behaviour  and  that  leads 
to  the  least  constraint  violation  possible.  The  inverse  design  variable  formulation  in¬ 
troduced  above  is  a  simple,  low-order  approximation  that  does  a  good  job  of  approxi¬ 
mating  the  constraint  functions,  especially  for  the  structural  constraints.  However,  it 
is  not  the  most  conservative  approximation  that  could  be  obtained  using  first-order 
expansions,  where  the  “most  conservative”  means  the  most  positive  for  constraints  of 
the  type  g(p)  <  0  and  the  most  negative  for  constraints  of  the  g(p )  >  0. 


The  hybrid  design  variable  approximation  obtains  the  most  conservative  possible 
approximation  by  combining  the  features  of  the  inverse  design  variable  approxima¬ 
tion  as  presented  above,  and  of  a  direct  design  variable  approximation,  where  the 
expansion  is  performed  with  respect  to  the  direct  design  variables  [Sta].  Consider 
approximating  the  constraint  function  using  the  first  two  terms  in  the  Taylor’s  series 
expansion  that  is  linear  in  the  direct  design  variables,  to  get 


<Up)  =  9(Po)  +  £ 


(Pj-Pjo) 


Po 


(3.8) 


or  by  an  expansion  that  is  linear  in  the  inverse  design  variables  t /_,,  as  defined  in  the 
previous  section,  as 


ff,(v)  =  ff(y0)+ 


(Vj  -  y10) 


Vo 


(3.9) 


Since  our  gradient  information  is  with  respect  to  the  direct  design  variables  pr  g t 


must  be  written  as 


v—'  2  dy> 
0 ,(p)  =  9(Po)  -  ojjr 

j=!  UI'J 


(3-1 0) 


A  criterion  for  determining  the  most  conservative  approximation  from  equation  (3.S) 
and  (3.10)  is  obtained  [Sta]  by  subtracting  equation  (3.10)  from  (3.8),  to  get 


f  1  dg 


(Pj-Pjof 


Po 


(3.11) 


S’nce  the  constraint  equations  used  in  the  problem  formulation  are  expressed  as 
g(x)  <  0,  it  can  be  seen,  with  reference  to  equation  (3.11),  that  g,  is  more  conservative 
(greater)  than  gL  when 


1  dg 
PjdPj 


P  0 


<  0 


(3.12) 


and  gL  is  more  conservative  (greater)  than  g,  when 


1  dg 
Pj  dp, 


'Po 


>  0 


(3.13) 


The  hybrid  approximation  formulation  is  motivated  by  the  desire  to  obtain  an 
approximation  that  best  predicts  the  actual  function  behaviour  without  violating  the 
constraint.  Then  let  the  hybrid  approximation  be  given  by 


9H(P)  =  9(P)  + 

j= i  °Pi 


Po 


(3.14) 


where 
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(?;  -  Pi.) 


^  = 


.,  1  dg  n 

if  — >0 
Pi  dPi  „ 

Fo 


~)(Pj-Pjo)  if  <0 

Pl  )  Pl  VPl  „ 

Fo 


(3.15) 


With  this  hybrid  approximation,  a  given  constraint  function  may  have  a  linear  ap¬ 
proximation  with  respect  to  one  variable,  and  an  inverse  approximation  with  respect 
to  another  variable,  hence  the  term  “hybrid”.  If  the  pj  are  physical  design  variables 
such  as  areas,  they  are  restricted  to  be  positive,  and  the  term  1/pj  can  be  discarded 
from  the  conditions  in  equations  (3.15). 


If  some  of  the  pj  are  non-physical  design  variables  which  are  unrestricted  in  sign, 
they  could  become  close  to  or  equal  to  zero  at  some  stage  in  the  design  process.  In 
such  a  situation,  an  inverse  design  variables  approximation  becomes  non-sensical,  and 
the  direct  design  variables  approximation  should  be  used.  Accordingly,  the  general 
hybrid  approximation  definition  for  Aj  in  equation  (3.15)  should  be  modified  to 


(Pj-Pio) 


*i  = 


.,  1  dg  i  i  » 

if  —  -r—  >0  or  \pj  <  A 

Pi  rj 

Fo 


“)(P;-Pjb)  if  <0 

Pj  )  Pj  “Pj  n 

Fo 


(3.16) 


where  A  is  a  “small’  number  specified  by  the  user. 

The  gradients  of  the  hybrid  approximation  function  can  be  found  by  differentiating 
equation  (3.14)  with  respect  to  the  direct  design  variables  p;,  to  get 


dgH{p )  _ 

dpj 


if  —  jr-  >0  or  |pj|  <  A 
Pi  dPi  jj 

F  o 


P*\  <h_  if  1  dg 
Pi  )  dPi  n  1  Pi  dPi 


(3.17) 
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3.5  Move-limits 


Move-limits  are  imposed  on  the  design  variables  during  each  approximate  problem 
solution.  This  is  done  in  an  attempt  to  restrain  the  design  variables  to  a  region  in 
which  the  explicit  function  approximations  remain  reasonably  accurate.  However, 
deciding  how  to  impose  these  move-limits  is  a  non-triviai  task.  The  local  curvature  of 
the  design  space  (i.e.  how  nonlinear  ar3  the  actual  constraint  surfaces  in  the  region 
about  the  expansion  point  of  the  approximations)  will  determine  the  move-limits,  with 
more  strict  move-limits  applied  in  regions  of  high  curvature,  and  less  strict  move-limits 
imposed  in  regions  of  low  curvature.  Since  second-derivative  information  is  required  to 
estimate  curvatures,  and  since  such  evaluations  are  very  expensive  computationally, 
imposing  move-limits  is  usually  reduced  to  an  art  based  on  past  experience.  It  is 
possible  to  obtain  the  second  derivatives  using  only  first  derivative  information,  and 
optimization  algorithms  that  do  this  are  termed  quasi-Newton  methods.  However, 
these  methods  typically  take  N  iterations  to  fill  the  Hessian,  and  can  be  very  costly 
if  N  is  large. 

For  the  purpose  of  this  work,  a  move-limits  factor  7  is  imposed  in  an  exponential 
form.  If  the  current  design  variable  and  approximation  expansion  vector  is  p0,  then 
the  upper  and  lower  bounds  on  the  design  variables  for  the  current  approximate 
problem  are  defined  as 


Pu  =  7Po 


(3.18) 


where  7  >  1.  The  limits  specified  in  equation  (3.18)  must  be  imposed  element  by 
element.  Note  that  since  the  design  variables  in  this  example  will  be  structural  design 
variables  only,  they  are  restricted  to  be  positive.  Obviously,  equation  (3.18)  must  be 
modified  if  the  design  variables  can  be  negative.  The  exponential  form  of  the  move- 
limit  factor  is  defined  by  the  particular  choice  of  7min  and  ymax>  as  shown  in  Figure  3.2. 
The  move  limits  get  more  restrictive  as  the  optimization  progresses,  and  as  the  design 
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supposedly  gets  closer  to  its  optimum  point.  Since  the  initial  design  may  be  very  far 
from  the  optimum  point,  allowing  large  moves  in  the  first  few  iterations  enables  the 
design  to  jump  around  the  design  space  somewhat  during  the  initial  stages  of  the 
design  process. 

Initially,  with  7  large,  the  design  is  more  likely  to  be  trapped  in  regions  of  a  deep 
local  minima  “well”,  and  is  more  likely  to  escape  from  shallow  local  minima  wells. 
As  the  move- limits  become  more  restrictive,  i.e.  as  7  — ►  1,  the  design  will  become 
trapped  in  a  local  minimum.  This  process  is  not  unlike  the  process  of  “simulated 
annealing”  [Kir-2],  although  obviously  not  as  general  or  controllable.  Unfortunately, 
there  are  no  general  conditions  which  will  indicate  which  local  minimum  is  the  global 
minimum,  or  where  it  can  be  found,  except  in  cases  where  the  design  space  is  convex. 
The  design  space  for  this  problem  is  definitely  not  convex. 

3.6  Scaling  to  the  Controller  Constraints  in  the 
case  of  Full  State  Feedback  Control 

In  the  case  when  full  state  feedback  is  to  be  used,  it  was  seen  in  Chapter  2  that  it 
is  reasonable  to  assume  that  the  controller  at  each  step  in  the  design  process  is  an 
appropriate  LQR  controller.  In  this  case  then,  it  becomes  practical  to  scale  the  struc¬ 
ture  to  the  closest  control  effort  constraint  and  closest  output  response  constraint 
simultaneously.  The  variables  with  which  the  structure  is  scaled  are  the  structural 
design  variables  (elemental  areas  or  thicknesses),  and  the  Lagrange  multipliers  asso¬ 
ciated  with  the  two  controller  constraints  Au  and  Ay  (where  for  clarity,  and  without 
loss  of  generality,  the  subscripts  on  the  A’s  that  refer  to  the  particular  control  effort 
or  output  response  constraints  under  consideration  have  been  dropped). 

Note  that  changing  the  values  of  Au  and  \y  cannot  independently  change  the  values 
of  uma  =  tr(I\TRKX)  and  ym,  =  tr(HjWHdX),  because  in  the  LQR  problem,  only 
the  ratio  of  Au  to  Av  is  important.  One  can  choose  the  ratio  (Au/Ay)  to  satisfy  one  of 
the  control  constraints  —  say  ums.  Then  yms  will  not  in  general  be  satisfied.  Suppose 
ym,  is  too  large  (i.e.  ym3  >  a2)  at  the  particular  point  where  ums  is  satisfied.  Then 
the  only  way  one  can  satisfy  the  yma  constraint  is  to  increase  the  sizes  of  at  least  some 
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of  the  structural  members. 


It  makes  sense  that  one  needs  to  change  the  structure  to  satisfy  both  control  con¬ 
straints.  If  the  control  constraints  could  be  satisfied  by  simply  choosing  appropriate 
controller  parameters,  then  there  would  be  no  interaction  between  structural  opti¬ 
mization  and  controller  optimization.  Intuitively,  it  can  be  seen  that  this  is  not  the 
case.  In  the  remainder  of  this  section,  a  method  is  developed  to  scale  both  the  ratio 
(Au/Ay),  and  the  cross-sectional  areas  of  truss  members,  to  exactly  satisfy  both  mean 
square  control  constraints. 

Note  that  each  member  of  the  structure  will  be  scaled  by  the  same  amount  to 
fulfill  our  goals.  Obviously,  this  method  is  not  absolutely  mandated,  and  some  other 
approach  could  be  used  where  the  design  variables  are  not  scaled  equally.  However, 
this  would  then  be  resizing  rather  than  scaling ,  a  process  normally  left  to  the  nonlinear 
programming  algorithm.  Assume  that,  at  iteration  i,  values  for  (uma),,  (ymj),,  the 
(Pj)»  f°r  j  =  1, . . . , (number  of  structural  members),  and  (Au),  (or  equivalently,  the 
ratio  (Au/Ay)i),  are  known.  Now  define 


(K)i 

(Att)i-i 

=  A, 

(3.19) 

( Pi)i 
(Pj)i-i 

=  Si 

(3.20) 

(«».»), 

=  Ui 

(3.21) 

(l/m*);-! 

=  Vi 

(3.22) 

Our  final  scaling  aim  is  to  set  umj  =  /?2  and  yma  =  a2.  Therefore,  let 


(3.24) 


(l/mj)t 

a2 

Further,  it  is  assumed  that  the  mean  square  values  (um,),  and  ( ym3)i  will  change,  as 
a  result  of  changes  to  (Au),-  and  the  (pj)i,  according  to  the  equations 


«.-+i 


Vi+ 1 


(^mj)i+l 

{Um3  )i 

=  A3’  8b> 

«+l  »+l 

(3.25) 

{Vma)i+ 1 
(2/ma)i 

=  ACl  8d' 

>  +  l  •  +  ! 

(3.26) 

where  <q,  b{,  C{  and  dt  are  constants. 


One  would  first  use  initial  (educated)  guesses  for  eq,  b{,  c,  and  d,  to  find  the 
appropriate  changes  to  (Au)0  and  the  (pj)o  required  to  satisfy  u i  =  /3 2  and  rq  =  a2. 
The  actual  changes  to  u0  and  y0  due  to  the  initial  values  of  a0,  b0,  cq  and  d0  will 
probably  not  be  those  predicted,  and  the  process  will  have  to  be  repeated  in  an 
iterative  fashion  to  obtain  the  desired  result.  However,  by  using  the  actual  changes 
to  iq  and  jq  at  each  step,  updated  and  improved  values  for  the  constants  a;,  6;,  c, 
and  d,  can  be  estimated  from  the  previous  calculations.  There  are  two  steps  to  each 
iteration. 


(a)  Evaluating  Al+1  and  <51+1  given  a,,  c<  and  d, 

We  desire,  at  the  next  step,  (umj),+1  =  02  and  (j/m,),+i  =  a2.  That  is, 


u.+i 


(^mi)i+l 


-l 

U 

•  +1./3 


(3.27) 


yi+ 1 


{llma)i+ 1 
(ymj)i 


=  Ac'  8d>  = 

>+i  >+i 


-i 

y 


(3.28) 

(3.29) 


These  provide  2  equations  in  the  2  unknowns  At+1  and  6  .  The  two  unknowns  can 
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be  found  by  taking  logarithms  of  equations  (3.27)  and  (3.28),  to  get 


log(A,+1 )  =  - 

di  l°g(ui+1  „)  —  Mog(yi+1  a) 

(3.30) 

a.idi  —  b{Ci 

losK+1)  =  ■ 

-Ci  log(ui+li0) 

(3.31) 

d\d\ 

which  will  be  valid  as  long  as 

atdi  -  biC,  ±  0  (3.32) 

If  this  condition  is  violated,  this  would  indicate  that  the  proposed  changes  cannot 
independently  scale  uma  and  yms.  However,  it  must  be  remembered  that  the  constants 
a*,  6,,  C{  and  di  used  here  are  not  exactly  known,  so  the  condition  (3.32)  may  sometimes 
be  violated  because  of  the  particular  values  of  the  constants  at  that  time.  Numerical 
experience  with  this  solution  algorithm  indicates  that  condition  (3.32)  is  not  restictive. 
If  this  condition  does  occur  at  some  point,  then  the  unknowns  are  set  to  their  values 
on  the  previous  iteration,  and  the  determinant  is  unlikely  to  be  zero  again  on  the 
next  iteration. 

(b)  Evaluating  at+l,  bi+l ,  c,+1  and  di+l  given  Ai+1,  5j+1,  u1+1  and  y,+l 

The  only  restriction  on  the  initial  values  a0,  b0,  Co  and  d0  is  that  as  defined  by 
equation  (3.32).  However,  one  can  obtain  initial  educated  guesses  by  performing 
parameter  tests  on  the  system  if  one  so  desires.  After  the  first  iteration,  there  is  only 
one  back  point  to  use,  hence  only  two  of  the  four  parameters  can  be  updated.  In  this 
work,  the  parameters  updated  are  a0  and  do,  while  setting  £>i  =  b0  and  ct  =  cq.  We 
know  A,,  <5,,  it,  and  y%,  and  wish  to  fit  a!  and  tfj  to 


=  A  “5<5f° 

(3.33) 

=  Af  6*' 

(3.34) 
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Again,  taking  logarithms  of  equations  (3.33)  and  (3.34)  gives 


logK )  ~  bo  togto  ) 
M  A, ) 


(3.35) 


_  log^J-colo^A,) 
1  log(^) 


(3.36) 


After  the  second  iteration  there  are  back  points  available,  so  all  of  the  parameters 
can  be  updated.  Consider 


y  =  Aa,+1«56,+1 

•+1  »+l  i+1 

(3.37) 

y  =  Ac,+1<5rfi+1 

y.+l  “i+l  Wi+I 

(3.38) 

u  =A  a*+,66t+> 

•  *  1 

(3.39) 

y  =  ACi+16*+1 

*'»  it 

(3.40) 

(3.41) 

and  solving  equations  (3.37)  through  (3.40)  for  the  unknowns  gives 


=  log(£, )  logK+, )  -  logto+t )  logK ) 
,+1  log(A,+1)log(5,)  -  log(AJ  log(£j+1 ) 

i°g(A,+1)iog(tt,)  ~  MAJiogK+J 

1+1  M Aj+1 )  l°g(5. )  -  log(- )  log(^+i ) 

log(S, )  log(y^, )  -  logto^ )  log(y, ) 

',+1  log( A,+1 )  log(«5, )  -  iog( A, )  log(£j+1 ) 

f  _  1<?g( A,4i )  log(y. )  -  los( A. )  Iog(y,.H ) 
,+1  M  Aj+1 )  log(6, )  -  log(  A, )  log(<5)+, ) 


(3.42) 

(3.43) 

(3.44) 

(3.45) 
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which  will  be  valid  as  long  as 

log(A„)log(«,)  -  tog(A,)log(6,„)  0  (3.46) 

Once  again,  if  this  condition  is  violated,  then  the  unknowns  should  be  set  to  their 
values  on  the  previous  iteration,  and  the  iteration  repeated. 

3.T  Structural  Analysis  and  State  Space  Formu¬ 
lation 

Large  flexible  structures  are  distributed  parameter  systems;  that  is,  their  mass,  damp¬ 
ing  and  stiffness  characteristics  are  described  by  variables  depending  on  time  and 
space.  Therefore,  the  governing  equations  of  motion  will  be  partial  differential  equa¬ 
tions,  which  are  theoretically  of  infinite  dimension.  To  model  these  structures  in 
a  way  which  facilitates  easy  and  efficient  solutions  to  the  equations  of  motion,  the 
structures  are  normally  discretized,  commonly  by  using  the  well  known  finite  element 
technique  (NASTRAN  for  example).  Discretization  reduces  the  structure  to  a  finite 
degree  of  freedom  system,  where  the  governing  equations  of  motion  will  be  ordinary 
differential  equations  of  finite  dimension. 

In  this  work,  a  finite  element  model  of  the  structure  is  formed,  and  only  structures 
which  are  a  collection  of  truss  elements  are  considered.  The  problem  formulation  and 
solution  techniques  employed  are  general  however,  and  the  extension  to  structures 
composed  of  beam  or  plate  elements  should  be  straightforward.  A  truss  element 
allows  only  axial  forces  and  displacements,  and  for  a  fixed  material  and  geometry,  the 
only  design  parameter  is  the  cross-  sectional  area.  Other  elements  may  exhibit  more 
design  freedom.  For  example,  a  beam  element  may  have  cross-sectional  area  and  area 
moments  of  inertia  as  design  parameters. 

The  equations  of  motion  for  the  discretized  structure,  subject  to  control  forces  in 
a  disturbance  environment,  can  be  written  in  the  form 

Mi  +  ci  +  K,Z  =  Du  +  GWl  w  (3.47) 
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where  £  is  a  7-vector  of  generalized  displacements  ( q  is  the  number  of  degrees  of 
freedom  in  the  structure),  u  is  an  nu-vector  of  applied  control  forces,  w  is  an  nw- 
vector  of  disturbance  forces,  M  is  the  system  mass  matrix,  K,  is  the  system  stiffness 
matrix,  C  is  a  damping  matrix,  D  is  a  matrix  describing  the  locations  and  orientations 
of  the  force  actuators,  and  Gm  is  a  matrix  describing  the  points  of  application  and 
orientation  of  the  disturbances. 

The  system  mass  and  stiffness  matrices  A 1  and  Ks  are  built  up  from  informa¬ 
tion  about  the  elemental  material  properties,  geometry,  location,  orientation,  and 
boundary  conditions.  Knowledge  of  M  and  Ka  enables  calculation  of  such  open-loop 
structural  information  as  element  stresses,  nodal  displacements,  and  open-loop  nat¬ 
ural  frequencies.  The  damping  matrix  C  can  be  formed  during  the  finite  element 
process  using  special  damping  elements.  If  this  option  is  not  desirable  or  available, 
the  damping  effect  can  be  included  in  a  variety  of  ways  (see  Section  3.7.3). 

The  mass  and  stiffness  matrices  will  be  functions  of  the  free  elemental  parame¬ 
ters.  Furthermore,  these  matrices  will  be  linear  functions  of  these  structural  design 
variables  for  structures  comprised  of  only  truss,  membrane  and/or  shear  elements. 
For  the  case  considered  in  this  work,  where  our  structures  are  composed  of  only 
truss  elements,  the  structural  design  variables  will  be  the  cross-sectional  area  of  each 
element. 

Once  the  finite  element  analysis  has  been  performed,  and  the  mass  and  stiffness 
(and  possibly  damping)  matrices  obtained,  the  state  space  system  can  be  formed. 
This  involves  writing  the  second  order  equations  of  motion  (3.47)  as  a  first-order 
state  space  system  of  equations,  in  which  a  choice  of  “state”  must  be  made.  There 
are  many  such  states  that  could  be  chosen,  but  two  choices  are  the  most  common. 


3.7.1  Physical  State  Space  Model 

A  state  space  system  based  on  physical  variables  can  be  formed  by  defining  the 
n  —  27-dimensional  state  vector  x  as 
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X  = 


(3.4S) 


.k. 

The  elements  of  this  x  are  the  actual  degrees  of  freedom  of  the  elements.  Then,  the 
state  space  matrices  F,  G,  and  Gw  of  the  realization  of  equations  (2.1)  can  be  formed 
as 
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(3.49) 


A  general  design  output  matrix  Hd  will  be  of  the  form 

Hd -[^  Hi]  (3.50) 

where  Hx  gives  the  contribution  to  this  design  goal  from  the  displacement  states,  and 
H2  gives  the  contribution  to  this  design  goal  from  the  velocity  states. 

3.7.2  Modal  State  Space  Model 

An  alternative  choice  for  state  space  representation  is  to  first  form  the  equations  of 
motion  (3.47)  in  modal  form  [Mei],  as 

rj  +  +  [ui]ij  =  <&T  Du  +  $rGWlw  (3.51) 

where  T]  is  the  vector  of  generalized  modal  displacements  such  that  £  =  $17,  [u>?]  is  a 
diagonal  matrix  of  natural  frequencies  squared,  and  $  is  the  matrix  whose  columns 
are  the  mode  shapes  corresponding  to  the  natural  modes  of  vibration,  which  are 
assumed  to  be  normalized  such  that  =  I. 

To  form  the  state  space  realization  in  terms  of  these  modal  variables,  define  the 
state  vector  x  as 
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X  = 


(3.52) 


The  elements  of  this  x  are  the  generalized  displacements  and  velocities  of  the  modal 
coordinates  77.  Then,  the  state  space  matrices  F,  G,  and  Gw  of  the  realization  of 
equations  (2.1)  can  be  formed  as 


'  0  1  1  f  0  ' 

F  =  ,  G  =  ,  Gv 

-[w?] 


(3.53) 


Note  that  for  a  general  C  generated  from  passive  damping  element  information  (if 
this  is  even  possible),  the  transformed  damping  matrix  will  generally  not  be 

diagonal. 


A  general  design  output  matrix  Hj.  will  be  given  by 


Hd  =  [  Hx9  H2<1> 


(3.54) 


3.T.3  Damping  Models 


Damping  of  the  structure  is  difficult  to  model  in  the  traditional  finite  element  sense. 
Normally,  damping  in  a  certain  form  is  assumed  after  the  formation  of  the  mass  and 
stiffness  matrices.  One  form  is  Rayleigh  damping  [Bat],  where  the  damping  matrix 
is  formed  as  a  linear  combination  of  the  mass  and  stiffness  matrices  a s 


C  —  <zq M  -)-  a\Ks 


(3.55) 


If  the  damping  ratios  of  some  or  all  modes  are  either  known  or  to  be  specified,  C  can 
be  defined  using  the  Caughey  series  [Bat] 


C  =  M  Y,  ak{\r1K}k 


(3.56) 
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where  p  <  n  damping  ratios  i  =  1  ,  ...,p}  are  specified,  and  the  coefficients 
{a*,  k  =  1, . . .  ,p}  are  calculated  from  the  p  simultaneous  equations 

C.  =  2  +  a2<^f  +  •  •  ■  +  ap_iw2p  for  i  =  1, . . .  ,p  (3.57) 

Note  that  with  p  =  2,  this  result  reduces  to  the  Rayleigh  damping  as  in  equa¬ 
tion  (3.55).  For  large  p  in  equations  (3.57),  the  computational  effort  involved  in 
finding  the  coefficients  a*,  becomes  quite  large.  Also,  assuming  C  must  be  differenti¬ 
ated  at  some  stage,  the  Rayleigh  damping  as  in  equation  (3.55)  or  the  proportional 
damping  as  in  equations  (3.61)  is  in  a  much  better  form  than  that  of  equation  (3.57). 
A  disadvantage  of  Rayleigh  damping  is  that  the  higher  modes  tend  to  be  more  highly 
damped  than  conventional  wisdom  dictates. 

A  damping  model  where  the  damping  ratio  in  each  mode  is  the  same,  and  where 
the  modal  damping  itself  is  proportional  to  the  modal  frequency,  can  be  formed  by 
defining  C  as 


C  =  2<A'1/2  (3.58) 

where  7v'/2  is  the  symmetric  square  root  of  K  such  that 

KU2Km  =  I<  (3.59) 

and  (  is  the  damping  ratio  of  the  modes.  Note  that  with  the  C  as  defined  by  equa¬ 
tion  (3.58),  the  transformed  damping  matrix  in  a  modal  state  space  representation 
will  be 


$>TC<I>  =  2  C$TA'1/3$ 

=  2CM 

=  [20,]  (3.60) 


41 


where  [u>,]  is  a  diagonal  matrix  of  modal  frequencies,  and  [2£u\]  is  a  diagonal  damping 
matrix.  This  method  of  forming  a  damping  matrix  is  also  probably  not  practical  due 
to  the  effort  involved.  The  diagonal  nature  of  [2(wj]  means  that  the  modal  equations 
of  motion  (3.51)  will  decouple  with  this  form  of  damping  into  equations  of  motion  for 
each  mode.  However,  each  mode  is  restricted  with  this  form  of  damping  to  have  the 
same  damping  ratio  and  can  be  formed  by  simply  specifying  $TC$  to  be  a  diagonal 
matrix  [2£;c<;t].  That  is,  assume  proportional  damping  of  the  form 

=  2(iu>,6i]  (3.61) 

where  £,  is  the  modal  damping  parameter,  and  6,-j  is  the  Kronecka  delta.  The  as¬ 
sumption  in  proportional  damping  is  that  the  total  damping  in  the  structure  is  the 
sum  of  the  individual  damping  in  each  mode. 

For  state  space  realizations  based  on  physical  variables  in  this  work,  Rayleigh 
damping  as  in  equation  (3.55)  was  used  to  model  the  structural  damping.  For  state 
space  realizations  based  on  modal  variables  in  this  work,  damping  of  the  form  of  equa¬ 
tion  (3.61)  was  used  to  model  the  damping.  The  particular  values  for  the  constants 
a0,  «i,  and  £,■  used  are  given  later. 

3.8  Gradient  Analysis 

For  the  numerical  optimization  procedure  to  be  practical,  especially  as  the  dimen¬ 
sionality  increases  to  realistic  structures,  it  is  essential  that  it  be  possible  to  evaluate 
the  first-order  sensitivities  of  the  complex  objective  and  constraint  equations  in  an 
efficient  manner.  The  following  sections  summarize  the  analytical  derivation  of  the 
appropriate  gradients,  where  it  will  be  seen  that  such  efficient  expressions  can  indeed 
be  found. 

3.8.1  Objective  Function  Gradient 

For  structures  formed  from  truss/membrane/plate  type  elements,  the  objective  func¬ 
tion  (the  weight)  is  a  linear  function  of  the  finite  element  thicknesses  and/or  cross 
sectional  areas.  That  is, 
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(3.62) 


N,  Ne 

J(p3)  =  Y,  PiLip'i  +  Y.rrij 
<=i  j= i 

where  N,  is  the  number  of  structural  design  variables  (number  of  finite  elements, 
dimension  of  p3),  p,,  Li  and  p,  are  the  density,  fixed  geometric  size  (length  of  a 
truss  element,  or  area  of  a  membrane  or  shear  element)  and  the  structural  design 
variable  (cross-sectional  area  for  a  truss  element,  or  thickness  of  a  membrane  or  shear 
element)  of  the  ith  finite  element,  N(  is  the  number  of  nonstructural  lumped  masses 
on  the  structure,  and  m,j  is  the  jth  nonstructural  lumped  mass.  The  gradient  of  the 
objective  function  is  then  easily  constructed  as 

dJ(p3) 

-^  =  PiLi  for  t  =  1, ...  ,1V,  (3.63) 

which  is  constant  for  all  designs.  Note  that  the  gradient  of  the  objective  function 
with  respect  to  any  controller  design  variable  will  be  zero,  since  in  this  work  the  mass 
of  the  controllers  was  not  considered. 

3.8.2  Controller  Constraint  Gradients  for  the  case  of  Full 
State  Feedback  Control 

The  controller  constraints  are  given  by 


Oca  — 


tz{KTRiKiX} 

ft 


—  1  <  0  for  i  =  1 , . . . ,  rip 


9oc, 


_  tr[HjW,Hd,X}  1<Q 


af 


for  i  — 


(3.64) 

(3.65) 


where  X  is  the  covariance  matrix,  found  as  the  solution  to  the  Lyapunov  equation 


FclX  +  XFl  +  GwXwGl  =  0 


(3.66) 


where 
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Fa  =  (F  -  GI<) 


(3.67) 


For  this  particular  case,  it  is  assumed  that  the  feedback  gain  matrix  I\  at  any 
point  is  given  by  the  solution  of  the  appropriate  LQR  problem  as 

K  =  R~lGT\x  (3.68) 

FTAr  +  AXF  -  \XGR~1GT\X  +  W  =  0  (3.69) 


where 


J  1  ^u*  1  R  l 

(3.70) 

£jr)  HlW* 

(3.71) 

This  assumption  is  based  on  the  analysis  that  showed  that  the  first-order  Kuhn- 
Tucker  optimality  conditions  for  the  controller  variables  were  satisfied  for  an  LQR 
controller  of  this  type.  However,  this  result  only  holds  at  the  optimum  point,  and  not 
necessarily  everywhere  within  the  design  space.  This  introduces  another  constraint 
into  the  problem,  namely  the  LQR  constraint,  similar  to  the  covariance  equation 
constraint.  Using  this  approach  then,  both  X  and  K  are  now  treated  as  variables 
dependent  on  the  particular  values  of  the  structural  design  variables  p\ 

Therefore,  if  scalar  pj  represents  a  structural  design  variable,  the  controller  con¬ 
straint  gradients  are  given  by 


d{jcc,  1 

— —  —  -r^tr 

dp,  Pi 


—^-Rih'i  +  I<TR>^A  X  +  KjRiKi 
dP]  dp , )  dP] 
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(3.72) 


dfjoc, 

dp: 


at 


rtr 


(^w^+Kw^yx+HrWiHj£ 


(3.73) 


With  K  defined  as  in  equation  (3.68),  its  derivative  with  respect  to  pj  becomes 


dPj 


d-fA,  + 

OP:  dp} 


(3.74) 


where  can  be  found  by  differentiating  equation  (3.69),  and  using  equation  (3.68), 
to  be 


nr^A,  <5A- 

c'  +  aPj 


Fci  =  - 


'5F  _  dG  ' 
A  dPi  > 


Ar  +  A* 


'dF_  _  0G  > 
A  0Pj  ) 


+ 


3W 

dpi 


(3.75) 


where 


—  =  y  ( +  hTw^Il 
Sp,  kw)\dp,  '  * +  *  '3r. 


(3.76) 


Note  that  only  the  right-hand  side  of  the  Lyapunov  eqaution  (3.75)  changes  for 
differing  pj.  Since  most  of  the  effort  involved  in  the  numerical  solution  of  Lyapunov 
equations  is  in  either  a  decomposition  phase  (for  direct  solvers)  or  an  inversion  phase 
(iterative  solvers)  operating  on  the  left-hand  side  of  the  equation,  one  decomposition 
or  set  of  inversions  will  be  enough  to  solve  equation  (3.75)  for  all  pj. 

The  partitions  of  required  for  each  of  the  control  effort  constraint  gradients 

9p,  &pj  u 

as  in  equation  (3.72)  can  be  found  easily  knowing  the  partitions  of  K  used  to  define 
the  original  constraint  functions  (3.64).  From  equations  (3.72),  the  gradient  of  the 
covariance  matrix  is  required,  and  can  be  found  by  differentiating  equation  (3.66)  to 
get 


Fa 


0X_ 

dPi 


dX  T 

+  ^+n- 


(3.77) 
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where 


= 


0Fd  OFjt  dGw 
ld^x+x-^-+  dPJ 


XwGl  +  GwXt 


dGl 

'  dP] 


dFcl  =  (W_dG  dK\ 
dPj  \dPj  dPj  V  dPj) 


(3.78) 

(3.79) 


Note  again  that  only  the  right-hand  side  of  equation  (3.77)  changes  for  differing  Pj. 
Further,  equation  (3.77)  is  the  adjoint  of  equation  (3.75),  which  means  that  only  one 
decomposition  or  set  of  inverse  calculations  (which  will  of  course  need  to  be  stored) 
of  the  left-hand  side  of  equation  (3.75)  is  required  to  solve  both  equations  (3.75)  and 
(3.77)  for  all  Pj.  This  property  significantly  improves  the  efficiency  of  evaluating  the 
gradients  for  this  particular  problem.  However,  the  Lyapunov  equation  (3.77)  must 
still  be  solved  for  every  pj.  To  avoid  this,  define  the  matrices  Vx  and  Qx  as  the 
symmetric  positive  definite  solutions  to  the  respective  Lyapunov  equations  [Can] 


FjtVi  +  ViFcl  +  Kj  RiKi  =  0  (3.80) 

Fd  Qi  +  QiFcl  +  HjW,Hdt  =  0  (3.81) 

so  that 

K'T RJu  =  -  (FjP,  +  ViFcl)  (3.82) 

=  -  (FjlQl  +  QiFcl)  (3.83) 


Substituting  equations  (3. 82)  and  (3.83)  into  equations  (3.72)  gives 


do. 


cct 


dPj 


=  ^tr 


d^LR,K,+KjR,d-^)x-{FlV,  +  V,F^ 
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(3.84) 

dpoc,  _  1  . 

dp:  ~  *r 


®  + ^w-w x  -  w  * + ca)  i. 


Since  tr[.Ai?]  =  tr[BA]  for  any  square  AB ,  the  above  equations  can  be  rewritten  as 


^Pcci 

dp. 


RiKi  +  KfRi 


X -Vi 


(3.S5) 


dffoc, 

dpj 


— tr 
a? 


'dHj 

LI  dp. 


and  with  reference  to  equation  (3.77)  becomes 

[&**'+ *ra>ist;)x+Pt‘H‘. 

(3.86) 

{^w‘H*+Hlw‘it)x+e‘ni. 

The  advantage  of  equations  (3.86)  over  equations  (3.72)  is  that  only  the  two 
Lyapunov  equations  (3. SO)  and  (3.81),  for  Vi  and  Qj  respectively,  need  be  solved 
to  find  the  constraint  gradients  for  every  Pj,  rather  than  the  N  Lyapunov  equation 
solutions  required  to  evaluate  the  covariance  equation  gradient  for  every  pj.  Of  course, 
Hj  must  still  be  evaluated  for  every  pj,  but  this  would  have  been  required  anyway  in 
equation  (3.77). 

3.8.3  State  Space  Matrix  Gradients 

The  state  space  mati'ices  are  given  by  equations  (3.49)  and  (3.50)  for  a  realization 
based  on  physical  variables,  or  equations  (3.53)  and  (3.54)  for  a  realization  based  on 


ddoc,  _  1  . 

dpj  Qi 


9 Pec, 
dpi 
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modal  variables.  For  each  of  thses  cases,  the  gradients  of  these  state  space  matrices 
with  respect  to  the  structural  design  variables  is  required.  Obviously,  the  gradients 
of  these  matrices  with  respect  to  the  controller  design  variables  will  be  zero. 

Let  p3  represent  a  structural  design  variable.  Then  for  a  realization  based  on 
physical  variables,  the  gradients  of  the  state  space  matrices  are  given  by 
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(3.89) 


where 


dM’1 

dpj 


i  dM  ! 
=  -M  — — M 
dpj 


(3.90) 


is  the  derivative  of  the  inverse  of  the  mass  matrix  M  with  respect  to  the  structural 
design  variable  p:.  The  derivative  of  M  can  be  found  very  easily  since  it  is  simply 
the  derivative  of  the  local  elemental  mass  matrix  expressed  in  the  global  coordinate 
system.  In  the  case  when  physical  variables  are  used  as  above,  the  design  output 
matrices  Hdx  will  not  be  functions  of  the  structural  design  variables,  so  that 

^  =  0  (3.91) 

dpj 
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For  a  realization  based  on  modal  variables,  the  gradients  of  the  state  space  matrices 
are  given  by 


3F 

dpj 
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where 


(3.92) 
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(3.94) 


Ax  =  diag 
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1  dPi  dPi  J 
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(3.95) 

(3.96) 


The  design  output  matrices  Hd,  will  now  be  functions  of  the  structured  design  variables 
because  they  will  be  functions  of  the  eigenvector  matrix  $.  Therefore,  the  derivative 
of  Hd,  is  given  by 


dHd, 

dp. 


3$ 


H, 


3$ 


dp i  "dpj 


(3.97) 


For  the  case  when  modal  variables  are  used  as  above,  the  derivatives  of  the  eigen¬ 
values  and  eigenvectors  will  be  required.  This  poses  no  problem  if  the  instantaneous 
structure  has  no  repeated  eigenvalues,  and  the  relevant  expressions  can  be  found,  for 
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example,  in  [Fox,  Nel].  However,  most  researchers  in  the  past  seemed  to  avoid  cases 
where  repeated  structural  frequencies  were  present.  Only  very  recently  has  this  issue 
been  addressed  [Mil-2,  Jua],  and  methods  to  evaluate  the  eigenvalue  and  eigenvector 
derivatives  for  the  case  where  repeated  eigenvalues  are  present  are  given  in  Appendix 
B. 


3.9  Design  Variable  Linking 

Even  though  the  control  design  variables  have  been  removed  as  active  design  variables 
in  this  case,  the  number  of  structural  design  variables  can  still  be  very  large.  An 
optimal  solution  where  all  structural  elements  have  been  treated  as  independent  would 
probably  define  all  elements  as  being  a  different  size.  In  practical  terms,  it  is  not 
desirable  to  have  a  structure  where  every  element  is  different,  because  this  can  cause 
production  and  spares  problems.  A  process  known  as  design  variable  linking  can  serve 
to  reduce  the  number  of  different  sized  structural  elements  that  will  be  needed.  The 
simplest  case  of  linking  is  to  assign  a  single  variable  to  a  group  of  elements,  so  that 
all  elements  in  that  group  will  have  the  same  variable  value.  Note  that  this  process 
in  effect  introduces  more  constraints  into  the  design  process,  which  can  affect  the 
optimal  solution  found.  However,  a  trade  off  between  weight  and  ease  of  production 
may  be  necessary. 

If  there  is  no  a  priori  knowledge  of  the  optimal  structure  should  be  somehow 
symmetric,  it  is  difficult  to  incorporate  design  variable  linking  into  the  solution  pro¬ 
cedure.  Perhaps  a  better  approach  is  to  consider  the  linking  only  after  initial  results 
for  unlinked  designs  have  been  evaluated.  Regardless,  design  variable  linking  should 
be  carefully  considered  at  some  stage  in  the  design  process  for  any  practical  design. 
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Piguie  3.1.  Algorithm  for  Solution  by  Mathematical  Programming  and  Approxima¬ 
tion  Techniques 


5 


Chapter  4 

Example:  The  DRAPER  I 
Structure 


4.1  The  DRAPER  I  Structure 

The  DRAPER  I  structure  [Str],  which  has  been  used  as  a  generic  flexible  spacecraft 
model,  is  a  tetrahedral  truss  attached  to  the  ground  by  three  right-angled  bipods,  as 
shown  in  Figure  4.1.  Although  attached  to  the  ground,  this  model  will  act  as  a  typical 
flexible  structure  pointing  subsystem  (e.g.  antenna,  radar,  optical)  attached  to  a  rigid 
core.  Any  motion  would  then  be  with  respect  to  this  rigid  core,  and  transmit  forces  to 
it.  Consequently,  this  model  has  no  rigid  body  degrees  of  freedom.  The  finite  element 
model  has  12  truss  elements,  since  the  joints  are  pinned  and  transmit  no  moments. 
There  are  four  nodes  that  are  free  to  move  in  all  directions,  so  the  model  contains  12 
degrees  of  freedom.  The  structural  design  variables  axe  the  cross-sectional  areas  of 
each  of  the  12  truss  elements. 

In  [Str],  the  structure  is  set  up  as  a  non-dimensional  model.  However,  for  our 
purposes  we  use  a  dimensional  model  with  material  parameters  of  p  =  0.1  lb/in3, 
E  =  Young’s  Modulus  =  20  kpsi.  The  dimensional  values  E  and  p  were  chosen 
to  give  initial  numerical  values  of  structural  frequencies  for  the  dimensional  model 
roughly  comparable  to  those  of  the  non-dimensional  model.  The  nodal  coordinates 
are  listed  in  Table  4.1,  giving  the  length  of  the  six  upper  bars  as  10  feet  and  the 
length  of  the  lower  six  bipod  bars  as  feet  in  this  case.  The  structural  members 
are  numbered  as  defined  by  the  finite  element  model  connectivity  data,  summarized 
in  Table  4.2.  This  dimensional  model  contains  no  non-structural  mass. 
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Elements  7  through  12,  the  three  right-angled  bipods,  take  on  the  duties  of  velocity 
sensors  and  colocated  force  actuators.  Only  one  output  response  constraint  is  defined 
( na  =  1),  with  the  design  output  vector  yd  representing  the  line-of-sight  error  of  the 
top  vertex  \(x,y)  displacements  of  vertex  1].  The  disturbances,  labelled  w\  and  u>2  in 
Figure  4.2,  are  assumed  to  be  independent,  zero  mean,  Gaussian  disturbances  with 
intensity  1.0. 

Since  there  are  12  degrees  of  freedom  in  the  model  for  this  structure,  the  state- 
space  model  will  be  24t/l  order.  There  will  be  6  inputs  and  6  outputs,  corresponding 
to  the  6  legs  of  the  structure.  The  damping  added  to  the  state  space  system  will 
depend  on  the  state  space  realization  used.  For  cases  where  a  realization  based  on 
physical  variables  is  used,  the  damping  matrix  C  is  formed  to  be 

C  =  0.1M -b  O.OOltf 

For  cases  where  a  realization  based  on  modal  variables  was  used,  the  damping  ratio 
of  each  mode  was  specified  to  be  0.1  %  of  the  modal  frequencies  during  the  formation 
of  the  state  matrices. 

Only  one  control  effort  and  one  output  response  constraint  are  considered  in  this 
example.  The  weighting  matrices  R  and  W  are  set  to  the  identity  matrices,  so  that 
equal  weighting  is  given  to  all  components  of  u  and  yd.  The  minimum  cross-sectional 
areas  for  all  elements  was  specified  as  0.1  in2.  For  this  problem,  no  static  structural 
constraints  were  specified,  the  intent  being  to  investigate  the  effect  of  the  closed-loop 
controller  constraints  on  the  structural  design  optimization. 

4.2  Numerical  Considerations 

Initially,  nominal  values  of  all  elements  are  specified  arbitrarily.  For  differing  val¬ 
ues  of  allowable  expected  output  response  (a2)  and  allowable  expected  control  effort 
(t 3 2),  this  initial  design  is  scaled  by  adjusting  the  cross-sectional  areas  and  Lagrange 
multiplier  ratio,  to  meet  the  controller  constraints.  Using  the  gradient  information, 
the  design  variables  and  the  feedback  gains  are  computed  using  the  approximate  for¬ 
mulations  to  further  reduce  the  system  mass.  If,  upon  solution  of  the  approximate 
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problem,  constraints  are  violated,  the  structure  is  scaled  to  the  constraint  surfaces 
using  the  scaling  scheme  discussed.  This  will  change  the  weight  from  that  value 
returned  as  the  optimal  approximate  problem  weight.  Therefore,  even  though  the 
optimal  approximate  problem  weight  will  always  be  less  than  or  equal  to  the  weight 
before  the  approximate  problem  solution,  after  the  scaling  the  weight  change  for  that 
global  iteration  may  sometimes  be  positive. 

If  the  approximation  functions  are  a  high  fidelity  representation  of  the  actual 
constraint  surfaces,  then  during  the  solution  of  the  optimization  problem,  we  would 
not  expect  the  solution  to  stray  far  from  these  actual  constraint  surfaces.  That  is,  after 
each  iteration,  the  actual  constraint  values  would  be  close  to  the  values  predicted  by 
the  approximate  functions  In  practice,  this  is  only  true  when  the  approximate  problem 
domain  is  restricted  to  a  small  region  about  the  expansion  point  \  la  the  move-limits. 

During  the  first  few  global  iterations  of  the  solution  algorithm,  the  move-limits  are 
on  the  order  of  7max,  which  is  nominally  set  at  1000.  This  allows  virtually  unlimited 
movement  in  the  design  space.  After  several  global  iterations,  the  move-limits  are 
reduced,  in  the  manner  described  in  Section  3.6,  to  the  order  of  7 mtn.  This  parameter 
can  significantly  affect  the  convergence  properties  of  the  solution  algorithm,  and  is 
nominally  set  at  1.25,  allowing  the  design  variables  to  change  by  +25%  and  -20%  in 
one  iteration. 

Two  problems  may  occur  when  7  «  7mm-  First,  if  the  approximation  functions 
predict  the  actual  constraint  behaviour  very  well,  this  is  evidence  that  the  constraint 
surface  is  highly  linear  over  large  region.  In  such  a  case,  convergence  may  be  very 
slow,  and  would  be  improved  by  increasing  7min.  The  reverse  may  occur  if  in  this 
limiting  case  the  approximate  functions  predict  the  actual  constraint  behaviour  very 
poorly.  This  is  evidence  of  high  local  curvature  of  the  constraint  surfaces,  and  a 
reduction  in  the  move  limits  by  reducing  7m,„  will  generally  improve  convergence.  In 
fact,  obtaining  ultimate  convergence  may  actually  require  such  a  reduction. 

When  obtaining  optimal  designs  for  differing  and  a 2  values,  it  is  sometimes 
necessary  to  alter  the  value  of  7min  to  take  advantage  of  local  conditions  near  an 
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optimum  point.  This  can  be  accomplished  by  either  manually  stopping  the  solution 
process  and  changing  7 m;n  before  restarting,  or  by  developing  some  adaptive  online 
strategy  to  accomplish  the  same  task.  Changing  the  objectives  of  the  optimization, 
by  altering  a2  and  0 2  for  example,  will  often  require  a  change  in  7min.  The  design 
process  seems  to  be  insensitive  to  7 max,  which  is  kept  constant  for  all  cases  reported 
in  this  work. 

If  a  state  space  model  based  on  modal  variables  is  to  be  used,  the  open-loop 
frequencies  and  mode  shapes  will  be  required.  These  can  be  calculated  in  a  post 
processor  to  the  finite  element  analysis  code  for  given  mass  and  stiffness  matrices.  In 
our  program,  the  eigenvalues  are  first  estimated  using  the  Sturm  sequence  property  of 
the  eigenvalues  [Bat],  by  successive  bisection  of  the  range  containing  the  eigenvalues 
of  interest.  The  eigenvectors  of  interest  are  then  found  using  a  forward  iteration  tech¬ 
nique  after  applying  a  shift  of  the  estimated  eigenvalues  as  found  above.  Starting  with 
an  initial  trial  eigenvector  populated  by  ones,  the  iteration  technique  should  converge 
to  the  eigenvector  corresponding  to  the  eigenvalue  for  which  the  shift  approximates. 
Care  must  be  exercised  in  cases  where  “identical”  eigenvalues  are  detected.  In  these 
cases,  Gramm-Sclunidt  orthogonalization  is  used  at  the  beginning  and  end  of  the 
iteration  process  to  deflate  the  initial  and  final  values  of  the  eigenvectors  of  other 
eigenvectors  associated  with  the  repeated  eigenvalue.  Finally,  after  all  eigenvectors 
of  interest  are  extracted,  the  eigenvalue  estimates  are  updated  using  the  Rayleigh 
quotient  of  the  appropriate  eigenvectors.  As  a  consequence,  the  eigenvalues  will  be 
accurate  to  second  order  in  the  error  present  in  the  eigenvectors  [Nob]. 

The  forward  iteration  method  of  finding  eigenvectors  finds  one  eigenpair  at  a  time, 
and  consequently  numerical  errors  creep  in  and  are  transmitted  and  amplified  as  more 
and  more  eigenpairs  are  extracted.  Since  here  we  will  either  desire  the  entire  eigen- 
structure  for  smaller  systems,  or  a  set  number  of  the  lowest  frequency  eigenpairs  from 
which  to  construct  a  low-order  design  model  for  large  systems,  a  subspace  iteration 
technique  to  extract  all  the  desired  eigenpairs  at  once  would  probably  be  more  ap¬ 
propriate.  Errors  in  the  eigenvectors  can  be  very  detrimental  to  the  accuracy  of  the 
derivatives.  The  integrity  of  the  derivatives  is  very  important  in  any  gradient  based 
solution  approach. 
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It  is  well  known  that  the  evaluation  of  eigenvectors  is  numerically  difficult  and 
subject  to  numerical  error,  especially  for  large  systems  and  for  systems  with  repeated 
eigenvalues.  As  was  shown  previously,  the  eigenvalues  and  eigenvectors,  and  their 
respective  derivatives,  are  required  for  the  evaluation  of  the  derivatives  of  the  state 
matrices  with  respect  to  the  design  variables  in  cases  where  a  modal  state  space  repre¬ 
sentation  is  used.  Appendix  B  presents  two  methods  for  evaluating  these  derivatives, 
for  systems  both  with  and  without  repeated  eigenvalues. 

There  is  first  a  problem  numerically  in  deciding  when  repeated  eigenvalues  are 
present.  The  number  of  repeated  eigenvalues  detected  depends  upon  a  specified  de¬ 
tection  cutoff  value.  Two  eigenvalues  are  deemed  identical  if  their  relative  difference 
is  less  than  some  value  e,  which  is  specified  by  the  user.  The  choice  of  £  is  a  non-trivial 
task  since  numerically  two  eigenvalues  are  never  identical,  and  since  we  can  generally 
not  tell  apriori  the  number  and  multiplicities  of  repeated  eigenvalues  present  in  a  par¬ 
ticular  configuration.  For  example,  for  the  DRAPER  I  structure  with  all  structural 
design  variables  the  same  (equal  area  elements),  and  with  e  set  to  5  x  10-4,  four  pairs 
of  repeated  eigenvalues  are  detected.  However,  if  £  is  set  to  1  x  10-4,  only  three  pairs 
of  repeated  eigenvalues  are  detected. 

If  a  repeated  eigenvalue  is  not  detected  where  in  reality  there  is  one,  by  setting 
£  too  small,  then  the  eigenvector  obtained  will  be  treated  as  unique  when  in  fact 
it  is  not.  Derivatives  of  eigenvectors  for  repeated  eigenvalues  only  exist  for  specfic 
eigenvectors  from  the  associated  subspace.  Therefore,  failure  to  detect  a  repeated 
eigenvalue  can  result  in  large  discontinuities  in  the  associated  eigenvector  derivatives. 
Additionally,  the  eigenvector  derivatives  for  nonrepeated  eigenvalues  are  inversely  pro¬ 
portional  to  the  difference  between  that  eigenvalue  and  all  others  (see  Method  2  given 
in  Appendix  B).  Therefore,  a  non-detected  repeated  eigenvalue  will  dominate  the  cal¬ 
culated  derivative,  and  any  errors  in  it  will  be  greatly  amplified.  On  the  other  hand, 
if  a  repeated  eigenvalue  is  detected  where  in  reality  there  is  none,  due  to  setting  e  too 
large,  an  extra  eigenvector  will  be  included  into  the  basis  for  the  subspace  of  eigenvec¬ 
tors  associated  with  that,  particular  eigenvalue.  This  would  increase  the  dimension  of 
the  associated  subspacc  by  adding  a  mutually  orthogonal  eigenvector.  The  associated 
unique  eigenvectors  from  this  subspace  that  have  continuous  derivatives  would  then 
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include  the  effect  of  an  eigenvector  not  associated  with  the  repeated  eigenvalue. 

The  calculation  of  eigenvectors  and  eigenvector  derivatives  for  systems  with  re¬ 
peated  roots  is  a  particularly  troublesome  and  difficult  task.  Even  though  the  ma¬ 
jority  of  the  time  the  system  will  probably  not  have  repeated  eigenvalues  through 
symmetry,  this  can  occur  for  example  with  symmetric  structures  subject  to  sym¬ 
metric  disturbances.  In  determining  the  derivatives  of  eigenvectors  associated  with 
repeated  eigenvalues,  it  is  assumed  that  the  eigenvectors  associated  with  a  particular 
eigenvalue  form  a  subspace,  and  that  any  linear  coi-bination  of  eigenvectors  is  also 
an  eigenvector  associated  with  that  eigenvalue.  Numerically  however,  this  is  not  ex¬ 
actly  true,  due  to  the  errors  involved  in  finding  these  eigenvectors.  As  a  consequence, 
even  though  the  expressions  given  in  Appendix  B  for  derivatives  of  eigenvectors  as¬ 
sociated  with  repeated  eigenvalues  are  analytically  correct,  our  experience  has  shown 
that  these  derivatives  are  very  sensitive  to  errors  in  the  evaluated  eigenvectors.  These 
errors  are  carried  into  the  gradients  of  the  controller  constraints,  and  hence  into  the 
approximate  problem  generation.  When  at  an  active  controller  constraint,  solution 
of  the  approximate  problem  should  move  the  design  approximately  along  the  con¬ 
straint  surfaces.  However,  at  a  design  for  which  repeated  eigenvalues  are  detected 
this  may  not  occur  because  the  errors  in  the  gradients  compromise  the  integrity  of 
the  approximate  problem  at  this  point. 

One  possible  way  to  avoid  these  problems  is  to  use  a  state  space  formulation  based 
on  physical  variables.  Here,  only  the  mass  and  stiffness  matrices  are  required,  and 
these  are  calculated  during  the  finite  element  analysis  anyway.  Such  a  solution  to 
this  problem  is  not  practical  for  “large”  systems  however,  where  a  model  based  on  a 
subset  of  the  modes  will  usually  be  used  to  reduce  the  computational  and  memory 
storage  needs. 

4.3  Results 

Runs  were  made  using  CSOPT  on  the  DRAPER  I  structure  using  initially  an  inverse 
design  variable  approximation  for  all  constraint  functions.  The  initial  structure  was 
defined  with  all  structural  design  variables  set  at  10  in2,  and  with  the  Lagrangian 
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multipliers  Au  and  Ay  set  at  1.0.  This  set  of  initial  conditions,  summarized  in  Table  4.3, 
will  be  termed  the  symmetric  set  of  initial  conditions,  for  they  specify  a  structure  with 
a  number  of  vibrational  modes  of  the  same  frequency  (repeated  eigenvalues).  A  range 
of  allowable  expected  output  response  (a2)  of  1  x  10-5  in2  to  1  x  10-4  in2  in  steps  of 
1  x  10-4  in2,  and  allowable  expected  control  effort  (02)  of  50  lb2  to  SO  lb2  in  steps  of 
10  lb2  were  used. 

Tables  4.4  and  4.6  summarize  the  resulting  minimum  weight  in  pounds  found  by 
CSOPT  for  each  of  these  cases.  The  values  in  Table  4.4  correspond  to  the  case  where 
a  state-space  realization  based  on  the  modal  displacements  and  velocities  was  used, 
hereafter  called  CASE  A.  The  values  in  Table  4.6  correspond  to  the  case  where  a 
state-space  realization  based  on  the  physical  nodal  displacements  and  velocities  was 
used,  hereafter  called  CASE  B.  The  specific  value  of  7 m;n  used  in  each  of  the  cases 
listed  in  Tables  4.4  and  4.6  axe  given  in  Tables  4.5  and  4.7  respectively. 

We  would  intuitively  expect  two  trends  in  the  data  displayed  in  Tables  4.4  and 
4.6,  as  well  as  expecting  the  two  tables  to  be  the  same.  The  optimum  weight  should 
decrease  as  the  allowable  control  effort  02  is  increased  at  constant  allowable  output 
response  a2  (left  to  right  across  the  table),  and  the  optimum  weight  should  decrease 
as  the  allowable  output  response  a2  is  increased  at  constant  allowable  control  effort  02 
(down  the  table).  With  reference  to  Table  4.4,  we  can  see  that  this  trend  is  observed 
in  a  macroscopic  sense  only,  there  being  several  examples  where  this  trend  is  not 
observed.  For  example,  considering  the  first  column  of  Table  4.4,  which  corresponds 
to  0 2  =  50  lb2  for  varying  a2,  we  see  only  two  exceptions  to  the  expected  trends, 
these  being  at  a2  values  of  6  x  10"5  and  9  x  10-5.  Similar  results  are  observed  in  all 
other  columns  and  rows  of  Table  4.4. 

The  results  using  the  physical  variable  realization  are  more  consistent  than  using 
the  modal  variables,  although  still  not  totally  uniform.  In  Table  4.6,  we  find  that 
there  are  only  two  locations  where  our  expected  trends  do  not  seem  to  hold.  In  the 
case  where  0 2  =  50,  the  optimum  weight  for  a2  =  8  x  10-5  of  1115.2  lbs  is  actually 
lower  than  might  otherwise  be  expected  considering  the  other  values  in  this  column. 
The  only  other  exception  is  for  0 2  =  80  and  a2  =  7  x  10-5.  One  might  then  draw 
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the  conclusion  that  using  the  physical  nodal  variables  as  a  basis  for  our  state-space 
model  results  in  a  more  consistent  set  of  solutions,  with  respect  to  the  optimum 
weight  found,  than  does  the  case  where  the  modal  variables  are  used.  We  cannot 
however  expand  this  conclusion  to  say  that  one  state-space  model  produces  better 
results  (lower  optimum  weights)  than  the  other.  In  fact,  on  this  count  the  two  results 
are  very  inconsistent.  For  example,  for  fl2  =  50,  the  optinr  m  weights  found  in  CASE 
A  are  consistently  lower  than  those  found  in  CASE  B,  whereas  for  /32  =  80,  it  is 
the  CASE  B  results  that  are  consistently  better  (or  at  least  comparable).  It  might 
be  pointed  out  that  the  results  in  CASE  B  were  consistently  easier  to  obtain,  there 
being  no  need  to  alter  the  nominal  value  of  7m,n  (see  Table  4.7  c.f.  Table  4.5),  and 
the  number  of  global  iterations  required  for  convergence  being  consistently  lower. 

Some  understanding  of  these  contradictory  results  can  be  found  by  considering 
Tables  4.8  and  4.9,  which  give,  for  the  two  cases,  the  optimal  element  areas  found 
for  ft2  =  50  and  for  varying  a2.  Also  given  in  these  tables  are  the  number  of  global 
iterations  required  for  convergence,  the  final  values  of  the  Lagrange  multipliers  (which 
then  defines  the  LQR  controller),  and  the  initial  value  of  the  structural  design  vari¬ 
ables  (all  the  same  for  the  symmetric  set  of  initial  conditions)  at  which  the  initial 
scaled  system  satisfies  the  constraints.  Immediately  apparent  from  Table  4.8  is  a 
number  of  seemingly  separate  regions  of  the  design  space  into  which  this  structure 
has  converged.  For  example,  the  final  designs  for  a2  =  5  x  10~5  and  a2  =  7  x  10~5 
seem  to  be  similar  in  relative  structure.  Here,  “similar”  refers  to  the  relative  sizing  of 
the  structural  members,  in  that  design  variables  that  are  “larger”  in  one  design  are 
“larger”  in  the  other.  Both  these  designs  are  however  distinctly  different  from  those 
for  a2  =  1  x  10-5  and  a2  =  3x  10-5,  which  themselves  are  similar. 

Considering  Table  4.9  corresponding  to  a  physical  state  space  representation,  we 
see  a  totally  d’fferent  phenomenom.  Here,  the  solution  algorithm  seems  to  converge 
to  the  same  region  of  the  design  space,  with  the  exception  of  one  case  (a2  =  8  x  10-5). 
This  corresponds  to  one  of  the  cases  for  which  the  optimum  weight  does  not  fit  into 
the  pattern  suggested  by  the  other  cases  for  the  same  allowable  control  effort,  as 
mentioned  above.  This  particular  solution  obviously  lies  in  a  different  region  of  the 
design  space  than  do  the  other  solutions  corresponding  to  other  o2  values.  Also  note 
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that  the  number  of  global  iterations  required  for  convergence  is  very  low  in  all  cases 
except  the  one  for  which  another  region  of  the  design  space  is  encounted. 

The  conclusion  seems  to  be  that  we  are  converging  into  different  regions  of  the 
design  space  with  our  solution  algorithm,  and  that  there  are  numerous  local  minima. 
Several  columns  of  Table  4.8  seem  to  define  their  own  region  of  the  design  space, 
being  dissimilar  to  any  other  column.  In  other  words,  our  design  space  seems  to  have 
multidimensional  corrugations  leading  to  multiple  local  minima.  The  solutions  will  lie 
somewhere  on  the  intersection  hyperplane  between  the  surface  of  constant  allowable 
output  response  and  the  surface  of  constant  allowable  control  effort. 

This  corrugated  nature  of  the  design  space  can  be  illustrated  by  considering  the 
solutions  obtained,  for  the  same  constraint  case,  when  starting  from  different  initial 
conditions.  For  the  case  of  /?2  =  75  and  a2  =  1  x  10-5,  Tables  4.10  and  4.11  summarize 
the  results  of  runs  made  when  modal  state  space  realizations  and  physical  state  space 
realizations  are  used  respectively,  when  only  the  initial  conditions  are  varied.  The 
different  initial  conditions  are  defined  by  setting  all  structural  elements  equal  except 
the  first  (element  1),  to  which  is  added  a  percentage  of  the  size  of  other  elements. 
Even  with  this  limited  variation  in  the  initial  conditions,  there  are  seemingly  many 
distinct  regions  in  the  design  space  into  which  the  system  may  converge.  A  picture  of 
the  constraint  surfaces  as  a  one-dimensional  slice  of  the  multidimensional  space  will 
emerge  if  these  optimal  structures  are  varied  into  each  other  in  a  linear  fashion,  and 
the  constraint  values  are  calculated  between  each  case.  That  is,  the  structural  design 
variables  and  Lagrange  multipliers  axe  changed  linearly  from  the  optimal  values  in 
one  case  to  those  in  another  case.  Then  the  constraint  surfaces  obtained  would  be 
those  seen  when  travelling  in  a  straight  line  between  each  successive  point. 

The  results  of  such  an  analysis  are  shown  in  Figures  4.3  and  4.4  for  the  cases 
corresponding  to  those  given  in  Tables  4.10  and  4.11  respectively.  As  expected,  the 
weight  varies  linearly  between  the  cases,  but  it  is  the  constraint  curves  that  are  much 
more  revealing.  For  example,  considering  Figure  4.3,  one  can  see  that  between  case  1 
and  case  2,  there  is  a  “ridge”  of  output  response  larger  than  the  maximum  allowable 
value.  Similarly,  the  control  effort  first  decreases,  then  also  increases  to  a  ridge  of 
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high  value.  This  corresponds  to  a  hump  in  the  constraint  surfaces  between  the  two 
points  in  the  design  space.  Assuming  that  we  would  see  such  behaviour  when  moving 
in  every  direction  away  from  case  1  and  case  2,  rather  than  just  in  a  direction  between 
the  two  as  shown  in  Figure  4.3,  then  the  design  points  corresponding  to  these  cases 
would  represent  local  minima.  In  this  situation,  the  design  can  become  “trapped” 
in  the  such  a  locally  convex  region,  causing  the  solution  algorithm  to  converge  to 
different  points. 

With  reference  to  the  same  Figure  4.3,  one  can  see  that  both  the  output  response 
and  control  efforts  are  virtually  constant  between  cases  2  and  3,  while  the  weight 
increases  slightly  from  2053.0  lb  to  2090.6  lb.  This  indicates  that  case  2  and  case  3 
actually  represent  the  same  optimal  solution,  with  the  difference  being  accounted  for 
in  the  variance  allowed  by  the  convergence  criteria  used.  The  direction  in  the  design 
space  represented  by  the  movement  from  case  2  to  case  3  would  lie  in  the  intersection 
hyperplane  of  the  surfaces  of  constant  control  effort  and  output  response  constraints, 
and  would  be  at  a  shallow  angle  to  the  linear  surface  of  constant  weight. 

Similar  observations  can  be  made  considering  Figure  4.4,  where  none  of  the  cases 
appear  to  be  representing  the  same  optimal  solution,  although  the  solutions  are  very 
similar.  Cases  1  and  2  in  Figure  4.4  are  only  separated  by  a  low  ridge  of  constraint 
values,  whereas  cases  4  and  5  are  separated  by  a  high  ridge.  Although  it  is  difficult  to 
visualize,  and  impossible  to  sketch,  the  actual  (in  this  case)  13-dimensional  constraint 
surfaces  and  their  12-dimensional  surface  of  intersection,  travelling  between  specific 
points  in  the  design  space  reduces  the  surfaces  to  manageable  quantities.  Figures  4.3 
and  4.4  graphically  illustrate  a  design  space  that  is  a  very  complicated  function  of 
the  design  variables,  in  which  multiple  local  minima  abound. 

There  are  some  other  tests  we  can  make  on  the  hypothesis  that  we  are  becoming 
trapped  in  local  minima.  If  the  design  is  actually  trapped  in  a  local  minimum,  the 
solution  should  stay  in  the  vicinity  of  that  minimum  if  the  problem  is  changed  only 
slightly.  That  is,  if  a  converged  solution  is  used  as  the  initial  conditions  for  an 
optimization  run  where  the  constraint  objectives  are  changed  by  a  “small”  amount, 
then  the  new  problem  should  converge  to  a  point  that  is  “close  to”  the  initial  point. 
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Table  4.12  represents  such  a  situation.  Here,  the  solution  wa s  first  obtained  for  the 
case  where  (3 2  =  50  and  a2  =  1  x  10-5,  and  where  a  modal  state  space  realization  and 
inverse  design  variable  approximations  were  used.  This  converged  solution  was  then 
used  as  the  initial  conditions  for  the  cases  01  =  50  and  a2  =  2  x  10-s,  and  02  =  60 
and  a2  =  1  x  10-5.  Moving  down  each  column,  and  across  the  top  row,  of  Table  4.12, 
the  converged  solution  from  the  previous  case  was  used  as  the  initial  condition  for 
the  new  problem. 

As  can  be  seen  from  Table  4.12,  the  two  expected  trends  in  the  data,  as  mentioned 
previously,  are  now  observed  without  exception.  The  optimal  solutions  for  the  first 
two  columns  of  Table  4.12,  corresponding  to  cases  where  02  =  50  and  (32  =  60,  are 
given  in  Tables  4.13  and  4.14  respectively.  All  the  solutions  now  appear  to  be  in  the 
same  local  region  of  the  design  space,  as  evidenced  by  the  relative  sizing  of  the  optimal 
structures.  For  example,  note  that  in  all  converged  designs,  structural  elements  9,  10, 
and  12  are  at  their  lower  gage  limit  of  0.1  in2,  and  that  the  first  structural  element 
is  the  largest  by  fax.  These  results  test  the  hypothesis  that  designs  are  converging  to 
local  minima,  and  indicate  that  the  local  optima  are  real. 

Convergence  to  local  minima  occurs  regardless  of  whether  physical  or  modal  state 
space  representations  are  used  to  model  the  structure.  The  final  converged  solution 
when  a  physical  state  space  realization  is  used  should  also  be  optimal  for  the  same 
constraint  objective  case  if  a  modal  state  space  realization  is  used  instead.  Such 
a  situation  is  documented  in  Table  4.15.  The  first  column  in  Table  4.15  gives  the 
final  converged  solution  for  the  case  when  /32  =  50  and  a2  =  2  x  10~5,  and  when  a 
physical  state  space  realization  and  inverse  design  variable  approximations  were  used 
(c.f.  column  2  of  Table  4.S).  When  this  solution  was  used  as  the  initial  condition 
for  an  optimization  run  for  the  same  case,  but  when  a  modal  state  space  realization 
was  used  rather  than  the  physical  realization,  the  results  in  the  second  column  of 
Table  4.15  were  obtained. 

Note  that  the  initial  values  for  the  control  effort  and  output  response  constraints 
for  the  case  when  a  modal  state  space  realization  was  used  (column  2  of  Table  4.15) 
were  not  exactly  50  and  2  x  10-5  respectively  as  expected.  Although  the  differences 
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were  small,  they  were  numerically  significant,  and  underscore  the  differences  between 
physical  and  modal  state  space  representations  of  the  same  structure.  The  differences 
come  about  due  to  the  numerical  errors  involved  in  eigenvalue  and  eigenvector  calcu¬ 
lations.  Because  of  this  small  difference,  the  structure  is  scaled  up  by  a  small  amount, 
and  the  optimal  solution  is  given  by  the  variable  values  in  the  second  column  of  Ta¬ 
ble  4.15.  If  this  structure  was  then  the  initial  condition  for  an  optimization  run  where 
a  physical  state  space  representation  was  again  used,  the  results  obtained  are  given 
in  the  third  column  of  Table  4.15.  Once  again,  there  is  a  small  discrepancy  between 
the  initial  constraint  values  and  those  expected,  for  the  reason  given  above.  Quite 
obviously  however,  all  three  designs  given  in  Table  4.15  are  the  same  local  optimum 
point  despite  the  small  differences.  Once  again,  this  indicates  that  the  apparent  local 
optima  do  actually  exist,  and  are  not  figments  of  a  numerical  imagination. 

The  only  difference  between  the  results  from  CASE  A  and  CASE  B  runs  is  that 
modal  and  physically  based  state-space  models  respectively  are  used.  The  structure 
is  initially  symmetric,  resulting  in  repeated  eigenvalues.  The  final  structures  of  CASE 
B  listed  in  Table  4.9  are  also  symmetric  with  repeated  eigenvalues  (except  that  for 
a2  =  8  x  10~5).  None  of  the  optimal  structures  of  CASE  A  listed  in  Table  4.8  are 
symmetric.  In  all  our  studies  to  date  it  seems  that  the  solution  algorithm  in  terms  of 
physical  state  variables  seems  to  be  more  predisposed  to  retain  the  symmetry  than 
in  the  cases  when  a  modal  state  space  is  used. 

To  illustrate  -,hat  is  happening,  consider  the  iteration  histories  for  the  same  run 
from  CASE  A  and  CASE  B,  in  particular  the  case  when  (3 2  =  50  and  a2  =  2  x  10-5. 
The  iteration  histories  are  contained  in  Tables  4.16  and  4.17  for  CASE  A  and  CASE 
B  respectively.  From  Table  4.16,  we  can  see  that  during  the  first  global  iteration, 
when  there  are  four  detected  repeated  eigenvalues,  and  unlimited  movement  in  the 
design  space  is  allowed,  the  inverse  design  variable  approximations  used  here  are  very 
poor  approximations  to  the  actual  constraint  behaviour.  However,  during  the  second 
global  iteration,  when  we  also  allow  virtually  unlimited  movement  in  the  design  space 
but  now  have  no  repeated  eigenvalues,  the  approximations  are  much  better. 

These  results  are  typical  of  all  runs  in  CASE  A.  The  combination  of  using  an 
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inverse  design  variable  approximation  and  a  modal  state  space  representation  seems  to 
result  in  very  poor  approximate  functions  when  repeated  eigenvalues  are  present.  We 
believe  the  reason  for  this  is  due  to  the  errors  present  in  the  derivatives  of  eigenvectors 
associated  with  repeated  eigenvalues,  the  causes  of  which  were  outlined  in  the  previous 
section.  However,  even  when  there  are  no  repeated  eigenvalues  present,  we  can  have 
poor  approximations  due  to  cumulative  errors  in  the  gradients  combined  with  an 
inverse  design  variable  approximation  formulation. 

Using  a  physical  nodal  state  space  representation  does  not  require  calculation  of 
eigenvalue  and  eigenvector  derivatives,  only  the  derivatives  of  the  mass,  stiffness  and 
damping  matrices.  These  are  very  easy  to  obtain  given  the  element  informations  and 
the  global  and  local  elemental  coordinate  system  definitions.  The  solution  algorithm 
here  seems  to  work  very  well  and  converge  quickly,  even  in  the  case  when  repeated 
eigenvalues  are  present,  as  shown  in  Table  4.17.  Again,  the  results  indicated  in 
Table  4.17  are  typical  for  all  runs  in  CASE  B. 

All  of  the  results  presented  so  far  axe  for  cases  where  inverse  design  variable 
approximations  are  used  for  the  constraint  functions.  However,  if  hybrid  design  vari¬ 
able  approximations  are  used,  similar  results  and  trends  to  those  already  noted  are 
observed.  The  optimum  weight  found  by  CSOPT  when  hybrid  design  variable  ap¬ 
proximations  and  both  modal  and  physical  state  space  representations  are  used  are 
given  in  Tables  4.18  and  4.20  respectively  (hereafter  termed  CASE  C  and  CASE  D 
respectively).  Tables  4.22  and  4.23  give  the  optimum  designs  for  02  =  50  and  varying 
a 2  for  CASE  C  and  CASE  D  respectively. 

From  Table  4.22  (CASE  C),  it  is  apparent  that  using  hybrid  design  variable  ap¬ 
proximations  and  modal  state  space  representations  has  led  to  optimal  designs  which 
are  much  more  similar  than  was  the  case  when  inverse  design  variable  approximations 
were  used  (compare  to  CASE  A,  Table  4.8).  That  is,  the  optimal  designs  listed  in 
Table  4.22  for  the  most  part  seem  to  lie  in  the  same  region  of  the  design  space,  since 
the  relative  sizes  of  the  elements  are  similar.  OF  course,  there  are  still  exceptions,  for 
example  when  a2  =  7  x  10-5  and  a2  =  9  x  10-5.  These  two  designs  are  the  two  that 
might  seem  anomolous  when  considering  the  0 2  =  50  column  of  Table  4.18.  Final  con- 
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verged  weights  however,  seem  to  be  much  larger  when  using  hybrid  approximations 
than  were  obtained  using  inverse  approximations.  This  can  be  explained  by  the  fact 
that  the  hybrid  approximations  give  a  better  approximation  of  the  constraint  surfaces 
than  do  the  inverse  approximations.  Therefore,  less  scaling  to  the  constraint  surfaces 
need  be  done,  and  the  design  tends  to  converge  quickly  to  a  local  minima  close  to 
the  specified  initial  conditions.  Even  shallow  local  minima  will  trap  the  design  easily 
because  the  approximations  are  predicting  the  actual  constraint  behaviour  very  well. 
Of  course,  if  the  initial  point  is  in  an  advantageous  point  in  the  design  space,  lower 
weight  solutions  may  be  obtained  using  hybrid  approximations  than  were  obtained 
using  inverse  approximations,  as  seen  in  some  cases  in  Table  4.23. 

If  a  physical  state  space  representation  is  used,  then  using  either  the  inverse  or  hy¬ 
brid  design  variable  approximations  seem  for  the  most  part  to  give  optimal  designs  for 
varying  a2  and  (3 2  that  lie  in  the  same  region  of  the  design  space.  In  fact,  comparing 
Tables  4.6  and  4.20,  many  of  the  optimum  weights  are  basically  identical  regardless 
of  whether  inverse  or  hybrid  design  variable  approximations  are  used,  even  though 
the  actual  optimum  design  variables  may  not  be  the  same  (compare  Tables  4.9  and 
4.23).  Once  again  however,  referring  to  Table  4.23  (CASE  D),  there  are  exceptions, 
notably  for  a2  =  2  x  10-5  and  a2  =  6  x  10-5. 

Our  conclusion  is  that  this  procedure  has  achieved  only  limited  success.  Solu¬ 
tions  can  be  obtained,  but  no  guarantee  can  be  given  that  these  solutions  will  be 
“good”,  and  no  algorithmic  way  exists  to  decide  how  to  go  about  getting  a  better 
(lower  weight)  solution.  There  does  seem  to  be  some  advantages  to  be  gained  in  not 
attempting  to  be  as  accurate  as  possible  in  approximating  the  constraint  functions. 
The  corrugated  nature  of  the  design  space  makes  the  problem  very  difficult,  and  these 
difficulties  can  only  get  worse  as  the  dimensionality  increases. 
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Table  4.1:  Nodal  Positions  for  the  DRAPER  I  Structure 
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Element 

Initial  value  (in2) 

1 

10.0 

2 

10.0 

3 

10.0 

4 

10.0 

5 

10.0 

6 

10.0 

7 

10.0 

8 

10.0 

9 

10.0 

10 

10.0 

11 

10.0 

12 

10.0 

=  1.0 

A,  *1.0 

Table  4.3:  Symmetric  Set  of  Initial  Conditions 
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Table  4.4:  Optimal  weight  using  a  modal  state-space  realization,  the  symmetric  set 
of  initial  conditions,  and  inverse  design  variable  approximations.  (CASE  A) 


Table  4.5:  7 mm  used  to  obtain  the  values  in  Table  4.4  (CASE  A) 
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Table  4.6:  Optimal  weight  using  a  physical  state-space  realization,  the  symmetric  set 
of  initial  conditions,  and  inverse  design  variable  approximations.  (CASE  B) 
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„2  /  w  l  n-5 ' 


Scaled  DV 


Final  DV 
1 
2 

3 

4 

5 

6 

7 

8 

9 

10 
11 

12 


iter,  for 
convergence 


90.437 


125.893 

23.336 

45.632 

6.389 

10.973 

9.800 

17.627 

17.719 

0.100 

0.100 

18.206 

0.100 


73.091 

53.545 

15.268 

10.704 

5.375 

6.209 

8.610 

0.124 

13.250 

9.332 

0.100 

0.100 


3 

4 

52.212 

EljWlMETigCTll 

70.557 

49.717 

43.868 

36.664 

21.137 

23.151 

14.833 

35.217 

21.510 

5.599 

6.318 

4.957 

12.207 

4.213 

9.637 

6.953 

5.365 

4.914 

10.775 

5.638 

9.310 

8.802 

0.150 

9.901 

0.100 

0.100 

0.100 

0.711 

0.313 

8.307 

13.746 

9.473 

8.352 

0.100 

9.266 

0.100 

mm 

1.4S51 

0.8632 

mu 

1.0 

1.0 

15 

32 

32 

a2  (xlO-5) 


Scaled  DV 


Final  DV 
1 
2 

3 

4 

5 

6 
7 
S 

9 

10 
11 
12 


Au 

A, 


iter,  for 
convergence 


6 


36.918 


25.936 

30.650 


2.501 

0.245 

9.107 

7.744 

8.5S3 


.7522 

1.0 


21 


7 

8 

9 

10 

34.180 

31.973 

29.979 

28.594 

32.194 

9.723 

16.757 

20.331 

21.255 

27.425 

32.495 

38.508 

20.736 

28.618 

16.565 

11.434 

5.141 

4.083 

8.79S 

2.041 

7.713 

6.737 

4.348 

2.369 

4.260 

3.815 

6.692 

2.934 

S.299 

0.100 

0.100 

0.152 

9.363 

0.100 

9.62S 

0.885 

0.100 

6.4S2 

1.332 

3.242 

S.183 

S.323 

0.100 

5.016 

7.379 

8.218 

7.737 

2.532 

0.100 

6.782 

10.099 

0.100 

0.8544 

0.7937 

0.9117 

1.7622 

1.0 

1.0 

1.0 

1.0 

18 

1 

22 

22 

49 

72 


a2 

(xlO-5) 

Hit 

1 

2 

3 

4 

mam 

Scaled  DV 

108.789 

76.924 

62.S07 

54.392 

48.650 

Final  DV 

1 

84.594 

61.454 

51.563 

45.089 

40.826 

2 

84.594 

61.454 

51.563 

45.089 

40.826 

3 

84.306 

61.524 

51.036 

45.095 

40.707 

4 

5.546 

5.008 

4.317 

4.207 

4.121 

5 

5.454 

5.001 

4.512 

4.338 

4.095 

6 

5.454 

5.001 

4.512 

4.338 

4.095 

7 

0.100 

0.100 

0.100 

0.100 

0.100 

8 

0.100 

0.100 

0.100 

0.100 

0.100 

9 

0.100 

0.100 

0.100 

0.100 

0.100 

10 

0.100 

0.100 

0.100 

0.100 

0.100 

11 

0.100 

0.100 

0.100 

0.100 

0.100 

12 

0.100 

0.100 

0.100 

0.100 

0.100 

’vC-vb 

■ 

2.S898 

IKIU 

KKjUjgj 

i  ■  EH 

1.0 

lEI 

iter,  for 

5 

8 

7 

6 

6 

convergence 

a 

iO 

1 

O 

X 

)  . 

6 

7 

8 

9 

10 

Scaled  DV 

44.411 

41.117 

3S.461 

36.261 

34.101 

Final  DV 

1 

37.634 

35.206 

35.938 

31.444 

30.105 

2 

37.634 

35.206 

36.090 

31.444 

30.105 

3 

37.519 

35.082 

7.280 

31.387 

30.066 

4 

3.425 

3.669 

2.885 

3.580 

3.369 

5 

4.177 

3.703 

2.866 

3.494 

3.349 

6 

4.177 

3.703 

3.456 

3.494 

3.349 

7 

0.100 

0.100 

7.681 

0.100 

0.100 

S 

0.100 

0.100 

0.359 

0.100 

0.100 

9 

0.100 

0.100 

0.384 

0.100 

0.100 

10 

0.100 

0.100 

7.000 

0.100 

0.100 

11 

0.100 

0.100 

0.100 

0.100 

0.100 

12 

0.100 

0.100 

0.100 

0.100 

0.100 

Au 

i  reran 

2.90SS 

Bigii 

Bpgai 

Ay 

IKM 

1.0 

B 

B§i“*£ 

iHI 

iter,  for 

5 

6 

23 

8 

7 

convergence 


Table  4.9:  Optimal  design  variables  for  (32  =  50  in  CASE  B 


73 


Case  1  Case  2  Case  3  Case  4  Case  5 


For  all  j  t^I,  the  initial  conditions  are: 

Case  1:  pi  =  pj 
Case  2:  Pi  =  Pj  +  3% 

Case  3:  pi  =  Pj  +  6% 

Case  4:  pi  =  pj  4-  10% 

Case  5:  Pi  =  Pj  +  50% 

For  all  cases,  (Au)0  =  1.0.  (Ay)0  =  1.0 

Table  4.10:  Optimal  values  when  /32  =  75  .md  a2  =  1  x  10-5  for  differing  initial 
conditions,  when  using  a  modal  state  space  model 
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Case  1 

Case  2 

Case  3 

Case  4 

Case  5 

Final  Wt. 

2006.6 

1952.6 

1875.5 

Final  DV 

1 

52.392 

47.085 

28.699 

32.257 

63.124 

2 

52.392 

48.486 

64.095 

64.893 

39.146 

3 

53.498 

59.061 

55.951 

46.706 

41.904 

4 

2.894 

3.565 

3.370 

5.534 

3.865 

5 

2.933 

2.486 

4.261 

4.591 

0.901 

6 

2.933 

1.621 

3.601 

4.992 

3.438 

7 

0.100 

0.556 

0.100 

0.237 

0.102 

8 

0.100 

0.100 

0.407 

0.100 

0.102 

9 

0.100 

0.100 

5.331 

0.100 

0.100 

10 

0.100 

0.100 

0.100 

0.100 

6.512 

11 

0.100 

0.100 

0.100 

0.100 

6.910 

12 

0.100 

0.506 

0.100 

9.956 

0.100 

Au 

2.8362 

2.8044 

2.3156 

2.1606 

2.3039 

For  all  ^  1,  the  initial  conditions  are: 

Case  1 :  pi  =  Pj 
Case  2:  px  =  pj  +  3% 

Case  3:  Pi  =  Pj  +  6% 

Case  4:  px  =  p:  4-  10% 

Case  5:  px  =  pj  +  50% 

For  all  cases,  (Au)0  =  1.0,  (Ay)o  =  1.0 

Table  4.11:  Optimal  values  when  /32  =  75  and  a2  =  1  x  10~5  for  differing  initial 
conditions,  when  using  a  physical  state  space  model 
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Figure  4.3:  Constraint  surfaces  between  cases  listed  in  Table  4.10 
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Control  effort 


Weight  (lbs) 


Case 


<D 

CO 

e 

o 

a, 

05 

<U 


3 

a 

-k-> 

o 


r-95 


^90 


H85 


-80 

-75 
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l65 


Figure  4.4:  Constraint  surfaces  between  cases  listed  in  Table  4.11 
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Control  effort 


a2  (xlO-5) 

0 

2 

50 

60 

70 

80 

1 

2847.1 

2466.0 

2166.6 

1969.2 

2 

2045.9 

if  82.1 

1579.7 

1443.2 

3 

1685.6 

1478.8 

1327.5 

1197.2 

4 

1470.9 

1298.0 

1162.0 

1052.0 

5 

1325  2 

1174.2 

1049.3 

952.1 

6 

1217.0 

1083.2 

965.9 

880.0 

7 

1134.1 

1012.4 

901.4 

823.8 

8 

1065.7 

955.3 

849.4 

776.5 

9 

1014.5 

907.9 

806.4 

735.3 

10 

961.3 

861.0 

770.0 

700.2 

Table  4.12:  Optimal  weight  using  a  modal  state-space  realization  and  inverse  design 
variable  approximations,  where  the  initial  condition  for  each  case  is  the  converged 
solution  from  the  previous  case. 
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79 


1 


a2  (xlO-5) 


5 


Final  DV 
1 
2 

3 

4 

5 

6 

7 

8 

9 

10 
11 
12 


iter,  for 
convergence 


114.695 

15.953 

39.232 

5.278 

9.266 

9.229 

14.382 

12.023 

0.100 

0.100 

15.158 

0.100 


.5486 

1.0 


20 


82.224 

11.936 

28.262 

4.230 

6.601 

6.831 

10.519 

8.285 

0.100 

0.100 

10.676 

0.100 


68.219 

9.905 

23.448 

3.510 

5.47S 

5.670 

8.747 

6.865 

0.100 

0.100 

8.848 

0.100 


.5948 

1.0 


2 


59.872 

8.695 

20.579 

3.082 

4.809 

4.978 

7.686 

6.016 

0.100 

0.100 

7.760 

0.100 


54.153 

7.866 

18.614 

2.789 

4.351 

4.505 

6.966 

5.421 

0.100 

0.100 

7.010 

0.100 


cr  (xlO-5) 


8 


Final  DV 
1 
2 

3 

4 

5 

6 

7 

8 

9 

10 
11 
12 


49.952 

7.258 

17.170 

2.575 

4.015 

4.157 

6.433 

4.982 

0.100 

0.100 

6.461 

0.100 


46.684 

6.784 

16.047 

2.40S 

3.754 

3.887 

6.018 

4.637 

0.100 

0.100 

6.034 

0.100 


44.047 

6.402 

15.141 

2.273 

3.543 

3.668 

5.683 

4.355 

0.100 

0.100 

5.691 

0.100 


41.861 

6.086 

14.390 

2.162 

3.369 

3.488 

5.406 

4.114 

0.100 

0.100 

5.406 

0.100 


39.029 

6.565 

13.696 

2.445 

3.122 

3.280 

4.564 

3.225 

0.100 

0.100 

4.685 

0.100 
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state  space  realization 


modal 

physical 

2.0575 

1.943 

50.3162 

49.876 

2077.3 

2100.9 

72.075 

70.662 

53.401 

52.344 

15.355 

15.203 

11.849 

11.718 

4.579 

4.S56 

6.921 

7.317 

10.486 

10.259 

0.152 

0.125 

16.378 

15.S47 

11.299 

9.6S7 

0.100 

0.100 

0.100 

0.100 

2100.9 

2067.8 

1.4140 

1.4510 

3 

2 

Initial  output 
response  (xlO-5) 


Initial  control 
effort 


Initial  weight 


Final  DV 
1 
2 

3 

4 

5 

6 

7 

8 

9 

10 
11 
12 


Final  weight 


iter,  for 


convergence 


73.091 

53.545 

15.268 

10.704 

5.375 

6.209 

8.610 

0.124 

13.250 

9.332 

0.100 

0.100 


2077.3 


1.4574 


Table  4.15:  Optimal  designs  for  /? 2  =  50  and  a2  =  2  x  10“5  when  modal  and  physical 
state  space  realizations  are  used. 
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Table  4.16:  Iteration  history  for  ft2  =  50  and  a2  =  2  x  10  5  for  CASE  A 
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iter. 

no. 

initial 

weight 

7 

after  approx,  prob.  solution 

scaling  for 
constraint 

satisfaction 

number  of 
repeated 
eigenvalues 

weight 

ElyJWvJ 

(xl0-s) 

ESI 

923.6 

— 

— 

533.7 

28.10 

7.692 

4 

7105.1 

1000 

6933.2 

46.97 

1.461 

0.8402 

4 

5S25.1 

368 

1644.6 

97.13 

0.9261 

1.695 

0 

27S8.9 

136 

2287.5 

51.95 

2.092 

1.080 

0 

2469.9 

51.0 

2399.1 

50.14 

2.000 

1.004 

0 

5 

2409.2 

19.5 

2397.2 

50.04 

2.001 

— 

3 

6 

2397.2 

7.98 

2399.9 

50.00 

1.999 

— 

3 

2399.7 

3.72 

2395.8 

50.01 

2.000 

— 

3 

2395.8 

2.16 

2395.4 

50.01 

2.000 

— 

3 

1 

2395.4 

— 

— 

50.01 

2.000 

— 

3 

Table  4.17:  Iteration  history  for  /?2  =  50  and  a2  =  2  x  10~5  for  CASE  B 


a2  (xlO-5) 


1 

2 

3 

4 

5 

6 


0 

2 

50 

60 

50S6.7 

5222.9 

3545. S 

3717.4 

2943.0 

2719.9 

2606.4 

2299.7 

2342.9 

2223.2 

218S.7 

1995.8 

1372.0 

1302.3 

183S.8 

1720.4 

1094.5 

1620.1 

1649.4 

1556.5 

Table  4.18:  Optimal  weight  using  a  modal  state-space  realization,  the  symmetric  set 
of  initial  conditions,  and  hybrid  design  variable  approximations.  (CASE  C) 


a 2  (xlO  5) 


1 

2 

3 

4 

5 

6 


Table  4.19:  jmtn  used  to  obtain  the  values  in  Table  4.18  (CASE  C) 
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Table  4.20:  Optimal  weight  using  a  physical  state-space  realization,  the  symmetric 
set  of  initial  conditions,  and  hybrid  design  variable  approximations.  (CASE  D) 


Table  4.21:  7mtn  used  to  obtain  the  values  in  Table  4.20  (CASE  D) 
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a2  (xlQ~5) 


1 

2 

3 

4 

n 

Scaled  DV 

90.437 

45.213 

40.441 

Final  DV 

1 

116.537 

81.263 

67.488 

59.206 

52.214 

2 

87.551 

63.283 

51.825 

45.141 

40.642 

3 

61.273 

42.350 

36.864 

31.630 

28.695 

4 

58.520 

39.845 

35.986 

30.204 

27.994 

5 

32.627 

23.613 

19.781 

15.285 

13.356 

6 

9.427 

6.661 

4.829 

4.572 

4.043 

7 

28.100 

Ini 

14.569 

15.321 

14.5S7 

8 

57.705 

Effil 

22.056 

25.805 

23.107 

9 

46.652 

32.261 

23.649 

24.185 

21.667 

10 

33.016 

20.119 

17.652 

21.179 

18.471 

11 

39.318 

25.292 

22.653 

23.688 

22.125 

12 

0.100 

0.100 

0.100 

0.101 

Au 

Eisa 

Kinm 

Hgn 

A, 

mm 

mi 

Hi 

1 

ms 

iter,  for 

30 

21 

18 

21 

17 

convergence 

a 

2  (xlO-5 

) 

■■ 

6 

7 

8 

9 

10 

Scaled  DV 

36.918 

29.979 

28.594 

Final  DV 

1 

47.435 

29.891 

41.464 

18.029 

37.093 

2 

37.598 

27.263 

31.927 

23.017 

2S.34S 

3 

26.607 

22.176 

22.521 

25.037 

20.347 

4 

28.581 

11.703 

21.748 

5.455 

19.389 

5 

10.641 

9.119 

10.791 

5.107 

9.535 

6 

3.741 

6.946 

3.188 

4.202 

2.903 

7 

12.179 

6.183 

10.791 

0.100 

9.717 

8 

24.414 

5.617 

18.756 

8.244 

16.051 

9 

21.866 

7.299 

15.194 

8.509 

15.714 

10 

18.713 

0.292 

14.S27 

0.100 

13.490 

11 

20.965 

6.080 

Sail 

9.849 

15.053 

12 

0.100 

0.100 

wm 

9.836 

0.100 

Au 

mm 

0.3781 

0.7831 

KjftnKl 

A, 

■El 

ms 

1.0 

iter,  for 

1  12 

33 

20 

21 

24 

convergence 

Table  4.22:  Optimal  design  variables  for  0 2  =  50  in  CASE  C 
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Q2  (xlQ-5) 


1 

2 

3 

4 

mm\ 

Scaled  DV 

108.789 

76.924 

62.807 

54.392 

48.650 

Final  DV 

1 

112.468 

25.786 

56.915 

49.314 

48.663 

2 

112.468 

28.096 

56.915 

49.314 

48.663 

3 

112.459 

101.308 

56.911 

49.310 

15.S32 

4 

0.842 

4.676 

13.774 

12.225 

3.667 

5 

0.842 

3.696 

13.784 

12.393 

4.634 

6 

0.842 

3.951 

13.784 

12.393 

4.634 

7 

0.100 

0.100 

0.608 

0.591 

1.125 

8 

0.100 

0.100 

0.604 

0.591 

0.100 

9 

0.100 

0.100 

0.604 

0.591 

0.100 

10 

0.100 

0.100 

0.608 

0.591 

1.123 

11 

0.100 

8.171 

0.598 

0.591 

0.100 

12 

0.100 

8.960 

0.598 

0.591 

0.100 

Au 

3.9224 

2.3377 

vmm 

2.1723 

EEEEU 

Av 

1.0 

1.0 

1.0 

1HS1 

iter,  for 

13 

32 

4 

4 

15 

convergence 

a‘ 

(xlO-5 

) 

6 

7 

8 

9 

10 

Scaled  DV 

44.411 

41.117 

38.461 

36.261 

34.401 

Final  DV 

1 

45.S80 

37.259 

34.749 

32.774 

31.123 

2 

45.880 

37.259 

34.749 

32.774 

31.123 

3 

12.703 

37.256 

34.745 

32.770 

Ena 

4 

4.152 

8.534 

8.733 

7.631 

7.326 

5 

3.771 

8.575 

8.736 

7.649 

7.344 

6 

3.766 

8.575 

8.736 

7.649 

7.349 

7 

0.970 

0.347 

0.400 

0.319 

0.313 

3 

0.100 

0.346 

0.383 

0.308 

0.303 

9 

0.100 

0.346 

0.383 

0.308 

0.303 

10 

0.970 

0.347 

0.400 

0.319 

0.313 

11 

0.100 

0.342 

0.393 

0.2S4 

0.283 

12 

0.100 

0  342 

0.393 

0.284 

0.2S3 

Au 

FgijjBjl 

2.3706 

gggj 

mmm 

A, 

1.0 

HIM 

WBm 

iter,  for 

16 

4 

6 

5 

5 

convergence 

Table  4.23:  Optimal  design  variables  for  /?2  =  50  in  CASE  D 
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Chapter  5 

Conclusions  and 
Recommendations 


The  problem  of  the  integrated  control/structure  design  optimization  of  large  flexible 
structures  has  been  investigated  for  the  case  of  full  state  feedback  control  of  structures 
subject  to  a  stochastic  disturbance  input  model.  Due  to  the  special  nature  of  the 
controller  and  problem  formulation,  the  controller  design  variables  (elements  of  K) 
were  found  to  be  given  by  the  solution  of  an  associated  LQR  problem.  Therefore,  these 
variables  are  not  treated  explicitly  as  design  variables  in  the  optimization  procedure, 
but  rather  sized  by  a  scaling  procedure  developed  to  identically  satisfy  the  constraints. 
A  research  computer  program  was  developed  to  solve  the  problem  using  a  sequential 
approximations  technique  and  commercially  available  nonlinear  programming  code. 
Very  efficient  expressions  for  the  required  gradients  were  derived  that  significantly 
reduced  the  computational  burden  of  the  solution  algorithm. 

The  solution  algorithm  was  applied  to  the  DRAPER  I  tetrahedral  truss  strucutre 
for  a  range  of  constraint  objectives.  Convergence  could  be  obtained,  and  an  optimal 
solution  found,  in  every  case  attempted  in  this  work.  However,  the  problem  was  found 
to  have  many  local  minima.  Although  this  is  generally  the  case  for  non-globally  convex 
problems,  it  is  the  dominant  factor  in  this  particular  problem.  The  solution  method 
was  very  sensitive  to  the  various  parameters  that  must  be  rather  arbitrarily  chosen, 
such  as  move-limits  and  initial  conditions.  Altering  either  of  these  two,  for  example, 
can  result  in  the  design  converging  to  distinctly  different  regions  of  the  design  space, 
with  a  resultant  difference  in  the  optimal  weight  found. 
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More  consistent  results  were  obtained,  when  varying  the  design  constraints,  for 
cruses  when  a  physically  based  state  space  realization  was  used  than  those  cases  where 
a  modal  state  space  realization  was  used,  although  often  the  optimal  weight  found 
was  larger.  Here,  “consistent”  refers  to  the  similarity  of  optimal  designs  when  the 
constraints  were  changed.  This  is  due  to  the  difficulties  encountered  in  evaluating 
accurate  derivatives  for  eigenvectors  associated  with  repeated  eigenvalues.  The  sen¬ 
sitivity  of  these  gradients  to  inaccuracies  in  the  eigenvectors  in  cases  where  repeated 
eigenvalues  are  detected  means  that  when  using  a  modally  based  state  space  system, 
the  solution  is  biased  away  from  symmetric  designs.  Cases  when  a  physically  based 
state  space  system  was  used  indicated  no  such  biasing,  and  in  fact  many  optimum 
designs  in  these  cases  exhibited  symmetry  (repeated  eigenvalues). 

There  was  a  great  deal  of  variation  in  the  optimal  designs  when  modal  state  space 
realizations  and  inverse  design  variable  approximations  were  used,  due  largely  to  the 
first  few  iteration  steps.  Initial  steps  tended  to  be  quite  large  due  to:  (i)  large  move- 
limits,  (ii)  the  inaccuracies  noted  in  the  repeated  eigenvector  derivatives,  and  (iii)  the 
significant  scaling  of  the  structure  due  to  inaccuracies  in  the  inverse  approximations. 
This  caused  the  design  to  “jump”  around  in  the  design  space  for  a  few  iterations,  not 
unlike  the  process  of  “simulated  annealing”,  and  the  optimum  weight  obtained  was 
more  often  than  not  lower  than  cases  where  the  initial  iterations  were  “well-behaved”, 
as  in  the  case  of  a  physical  state  space  representation  and/or  hybrid  design  variable 
approximations. 

For  this  reason,  and  because  of  the  apparent  existance  of  so  many  local  minima,  it 
is  recommended  that  future  work  should  investigate  various  options  for  resolving  the 
best  solutions.  Systematic  solutions  using  random  initial  conditions,  or  the  intentional 
inclusion  of  a  random  perturbing  factor  as  in  the  “simulated  annealing”  procedure, 
may  be  required  to  obtain  solutions. 

This  work  concentrated  on  the  case  when  full  state  feedback  control  is  used.  How¬ 
ever,  in  most  designs  of  practical  significance,  full  state  feedback  is  probably  not 
practical.  Future  work  should  include  studies  on  different  controller  methodologies, 
and  how  they  affect  the  solution  algorithm.  Examples  of  large  dimension  must  be 
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attempted  and  solved  to  give  confidence  that  this  solution  process  is  viable  in  realistic 
design  situations. 

Even  with  these  problems,  our  conclusion  is  that  the  combined  control/structure 
optimization  problem  has  a  high  payoff  potential,  and  should  continue  to  be  explored. 
The  ability  of  the  procedure  to  trade-off  structural  mass  verses  controller  energy  is 
an  important  tool  for  coupling  structural  system  and  controller  design  synthesis,  and 
will  be  crucial  in  the  design  of  realistic  structures  for  widespread  application  in  the 
coming  decades. 
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Appendix  A 


Results  Concerning  Matrix  Trace  Derivatives 


For  scalar  s  and  matrix  A  with  elements  a^,  define  the  matrix 


ds 

dA 


a s 


ds  _  ds 
dA  tj  “  daij 


(A.l) 


Then  the  following  results  on  the  trace  operator  hold: 


=  t 

(A.2) 

Atr!  AB]-±«{BA\  =  BT 

(A.3) 

A  tr[A'B]  =  ±  tr^.4]  =  B 

(A.4) 

A  tr[ZJT.4C]  =  BCT 

(A.5) 

A  tr[4rfl.-tl  =  (B+BT)/i 

(A. 6) 

A  tr[4B<4T]  =  A{B+Br) 

(A.7) 

At^,  =  2/1 

(A.S) 
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Appendix  B 

Derivatives  of  Eigenvectors  Associated  with  Repeated 
Eigenvalues  of  Structural  Systems 

This  appendix  presents  two  methods  for  obtaining  eigenvector  derivatives  of  struc¬ 
tural  systems  with  distinct  eigenvalue  derivatives.  The  first  method,  due  to  Mills- 
Currcn  [Mil-2],  requires  only  one  eigenvalue  and  its  associated  complete  subspace 
of  eigenvectors,  and  only  the  results  are  given  here.  The  second  method  is  derived 
here  and  assumes  knowledge  of  the  entire  system  eigenstructure.  For  a  more  general 
method  that  does  not  require  distinct  eigenvalue  derivatives,  the  interested  reader  is 
referred  to  Juang,  et.al.  [Jua]. 

Method  1 

The  real,  symmetric  structural  eigenprobiem  (defining  K  and  M  as  the  stiffness 
and  mass  matrices  respectively)  is 

Fi<f)t  =  0  for  i  =  l,..., n  (B.l) 

where 

Fi  =  [I\  -  A.M]  (B.2) 

with  mass  normalization 

<pj M<$>}  =  Sij  for  i  =  l,..., n  (B.3) 

where  8tJ  is  the  Kronecker  delta  function. 


98 


Let  the  first  and  second  derivatives  with  respect  to  some  parameter  be  denoted 
by  '  and  "  respectively.  Then  differentiating  equation  (B.l)  gives, 


F[<i}l  +  F,4>',  =  0  for  i  =  1 , . . . ,  n 


(B.l) 


Premultiplying  by  <f>J ,  and  using  the  orthonormality  condition  (B.3),  we  can  obtain 
an  expression  for  the  eigenvalue  derivatives  as 

A'  =  (f)J [A"  -  A.A/']^  for  i  =  1 , . . . ,  n  (B.5) 

Because  the  original  system  is  symmetric,  an  independent  eigenvector  will  exist 
for  every  eigenvalue.  However,  eigenvectors  associated  with  repeated  eigenvalues  will 
not  be  unique,  but  will  form  a  subspace  of  the  same  dimension  as  the  eigenvalue 
multiplicity,  from  which  any  vector  will  be  an  eigenvector  of  that  eigenvalue.  Equa¬ 
tion  (B.5)  cannot  be  used  to  find  the  derivatives  of  eigenvalues  with  multiplicities 
greater  than  one  until  the  arbitrary  nature  of  the  associated  eigenvectors  is  removed, 
as  the  derivatives  may  not  exist  for  every  choice  of  A. 

To  remove  the  arbitrariness  in  the  ip:,  an  eigenvector  continuity  condition  must 
be  imposed  to  insure  that  the  system  eigenvectors  remain  continuous  as  the  system 
is  perturbed  about  the  condition  where  the  repeated  eigenvalues  are  found.  This 
essentially  chooses  a  unique  set  of  eigenvectors  from  the  subspaces  associated  with 
the  repeated  eigenvalues,  any  ensures  existance  of  the  eigenvector  derivatives. 

Since  the  eigenvectors  associated  with  repeated  eigenvalues  form  a  subspace,  we 
may  write 


or 


=  Y, 

u  n. 


(B.6) 


(B.7) 
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In  the  above  equations,  the  4>t  are  the  eigenvectors  of  the  repeated  eigenvalue  A,, 
where  fe  Q,  =  {i,z  +  1, . . .  ,*  +  m,-  —  1}  and  m,  is  the  multiplicity  of  A,.  The  vectors 
rpj  are  linear  combinations  of  the  <fr(,  with  the  coefficients  a tj  to  be  determined.  The 
columns  of  contain  the  vectors  the  columns  of  contain  the  vectors  4>y  and 
a,j  is  the  ijth  element  of  A.  The  matrix  <J>,  is  the  matrix  of  nominal  eigenvectors 
associated  with  A,,  and  is  the  matrix  of  unique  eigenvectors  associated  with  A,. 
Substitution  of  equation  (B.7)  into  equation  (B.l)  shows  that  the  xj;J  satisfy  the 
original  eigenproblem  for  any  choice  of  the  a(y 

After  imposing  the  continuity  condition,  the  columns  of  A  and  the  respective 
eigenvalue  derivatives  are  found  by  solving  the  subeigenproblem 

[A"  -  A.M'jaj  =  A'a_,  (B.S) 

A  unique  matrix  A  can  be  found  in  this  manner  only  if  the  derivatives  of  repeated 
eigenvalues  are  distinct,  i.e.  iff 


a;  ±  a;  when  A,  =  A;  (B.9) 

Assuming  that  this  is  so,  the  A  matrix  then  defines  a  set  of  unique  linearly  inde¬ 
pendent  orthogonal  eigenvectors  for  the  original  eigenproblem  (B.l).  The  derivatives 
of  these  eigenvectors  can  now  be  found.  From  this  point,  assume  that  the  </>,  repre¬ 
sent  the  set  of  unique  eigenvectors  as  found  above.  Then  the  eigenvector  sensitivity 
is  assumed  to  have  the  form 


4>\  =  Vi  +  <J>,c,  for  i  =  1, . . . ,  n 


(B.10) 


where  V,  is  the  vector  solution  to 


FxVt  =  - f:4> , 


for  t  =  1, . . .  ,n 


(B.l  1 ) 


Since  A,  is  not  of  full  rank,  equation  (B.li)  cannot  be  solved  by  direct  inversion. 
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However.  :f  the  appropriate  rows  and  columns  are  eliminated  from  F,  along  with 
the  corresponding  rows  from  the  right-hand  side  of  equation  (B.ll),  then  a  valid 
solution  for  V  can  be  found.  Guidance  for  the  partitioning  of  equation  (B.ll)  can  be 
found  in  the  original  reference  [Mil-2],  and  will  not  be  reproduced  here. 

The  elements  c,j  of  the  vectors  Cj  can  be  found  from 


ca  =  +  2MV{)  for  i  =  l,..., n  (B.12) 


<t>] {!<"  -  2XiM'  -  \iM")<j)i  +  24>] FjVi 
2(A<  -  A') 


for  i,j  =  1,. 


n  and  i 


7^; 

(B.13) 


Method  2 

In  this  method,  we  again  assume  that  a  unique  set  of  eigenvectors  for  the  eigen- 
problem  have  been  found  as  outlined  in  Method  1  above.  However,  we  now  assume 
the  eigenvector  derivative  is  a  linear  combination  of  these  eigenvectors,  in  the  form 

#  =  =  (B14> 

3=1 


( 

► 


Substituting  equation  (B.14)  into  the  mass  orthonormality  condition  (B.3),  we  obtain 


(B.15) 


For  any  j  ^  i  such  that  A j  ^  A*,  expressions  for  the  coefficients  6|;  are  well  known 
[Fox],  and  are 


where 


(A,  -  A,) 


(B.16) 
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Z,  =  [!<'  -  A,  A/'] 


(n.  it) 


To  obtain  the  remaining  unknown  coefficients,  differentiate  equation  (B.4)  again, 
to  get 


F?<t>i  +  2 Ft#  +  Fitf  =  0  for  i  =  1, . . . ,  n 


(B.18) 


Premultiplying  equation  (B.18)  by  <f>J,  for  j  such  that  A,-  =  A j,  we  get 


<t>?F!'<t>i  +  2<t>jFl<f>'i  =  0 


(B-19) 


since  (j>J F,  =  0.  Substituting  equation  (B.14)  for  the  eigenvector  derivative,  we  get 


+  =0 


(B.20) 


Consider  the  first  term  of  equation  (B.20): 


=  <£j[A"'  -  2A'A/'  -  A, AT  -  A''A/]<£, 


=  [I<"  -2A'A/'- A,M"]  -  \"Sij 


(B.21) 


Now,  consider  the  second  term  of  equation  (B.20): 


20jF'(£M>,)  =  24^\K"  ~  KM  - 


=  24>t!{K"  ~  X\M  -  AjA/'] 
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2#(A"'  -  A'M  -  A, Ml  (£>,&)  + 

+  2(A'-A')<tfA/(£w,) 

2«jF'('f;i,(i0,')+26jj(A'-A;) 


(B.22) 


^=1 


To  reduce  the  first  term  on  the  right  hand  side  of  equation  (B.22),  substitute  equa¬ 
tion  (B.14)  into  the  transpose  of  equation  (B.4),  to  get 


(B.23) 


Since  \j  =  A,,  then  F}  =  F„  and 


0 

(A*  -  A 


if  A k  —  A, 
if  A*  ^  A, 


(B.24) 


Therefore, 


=  -2  t,  MV  -  V WfM  (  £  W  J 

tm\  \m=l  / 

<<n, 

=  -2  £  -  Ai)  (B.25) 


where  ft,-  =  {  k  :  A*  =  A,  }.  Since  we  already  have  expressions  for  b(j  and  ba  for 
t  ^ft,  given  in  equation  (B.16),  then  equation  (B.25)  becomes 


-2 


E(v-v> 


tml 


(tf  Zj»i)  (tfz.w 

(A,  -  A,)  (A,  -  A,)  . 
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(B.26) 


=  2£ 
/=1 
^n, 


(A,  -  A,) 


Substituting  equation  (B.26)  into  equation  (B.22),  then  equations  (B.21)  and 
(B.22)  into  equation  (B.20),  we  obtain 

v 


4>J[K"  -  2 A' A/'  -  -  A”<„  +  26,, (A'  -  A')  + 


+  2E 

tzt  1 


(A,  -  A,) 


(B.27) 


* 


Forcing  i  ^  j,  and  rearranging,  gives  the  required  coefficient  for  i,j  /  fl,  as 


£ 

/=1 

/f'n, 


(A,  -  A,) 


(B.2S) 


Note  that  differentiating  the  M-orthogonality  condition  (B.3)  between  <£,  and  (f>J 
gives 


tJ'Mti  +  4>]M4>'X  =  -<#  M'4>i  (B.29) 

Substitution  of  equation  (B.14)  into  equation  (B.29)  gives 

bij  +  bji  = —<f>J M'<f> ,  for  i  =  n  (B.30) 

Equation  (B.2S)  satisfies  equation  (B.30),  as  would  be  expected.  Also  note  that 
forcing  i  =  j  in  equation  (B.27)  gives  the  second  derivative  of  the  eigenvalues  as 
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K  =  2  £  + ^l/r  -  2A'iM'  - 


<•1 


f 

* 


(B.31) 


105 


®U.S.Oovarnm«nt  Printing  Office  1990— 748*056/24469 


